module shr_flux_mod

  ! atm/ocn/flux calculations

  use shr_kind_mod, only: r8 => SHR_KIND_r8, in => SHR_KIND_IN ! shared kinds
  use shr_const_mod
  use shr_sys_mod

  implicit none

  private


  public flux_atmOcn                ! computes atm/ocn fluxes
  public flux_atmOcn_diurnal        ! computes atm/ocn fluxes with diurnal cycle
  public flux_atmOcn_UA             ! computes atm/ocn fluxes using University of Ariz algorithm (Zeng et al., 1998)
  public flux_MOstability           ! boundary layer stability scales/functions
  public shr_flux_adjust_constants  ! adjust constant values used in flux calculations. (used by CAM as well)

  private psi_ua
  private qsat_ua
  private rough_ua
  private cuberoot
  private cor30a
  private psiuo
  private psit_30

  integer(in), parameter, public :: shr_flux_MOwScales   = 1 ! w scales  option
  integer(in), parameter, public :: shr_flux_MOfunctions = 2 ! functions option
  real   (r8), parameter, public :: shr_flux_MOgammaM    = 3.59_r8
  real   (r8), parameter, public :: shr_flux_MOgammaS    = 7.86_r8

  !--- rename kinds for local readability only ---

  integer, parameter :: debug = 0 ! internal debug level

  ! The follow variables are not declared as parameters so that they can be
  ! adjusted to support aquaplanet and potentially other simple model modes.
  ! The flux_adjust_constants subroutine is called to set the desired
  ! values.  The default values are from shr_const_mod.  Currently they are
  ! only used by the flux_atmocn routine.
  real   (r8) :: loc_zvir   = shr_const_zvir
  real   (r8) :: loc_cpdair = shr_const_cpdair
  real   (r8) :: loc_cpvir  = shr_const_cpvir
  real   (r8) :: loc_karman = shr_const_karman
  real   (r8) :: loc_g      = shr_const_g
  real   (r8) :: loc_latvap = shr_const_latvap
  real   (r8) :: loc_latice = shr_const_latice
  real   (r8) :: loc_stebol = shr_const_stebol
  real   (r8) :: loc_tkfrz  = shr_const_tkfrz

  ! These control convergence of the iterative flux calculation
  ! (For Large and Pond scheme only; not UA or COARE).
  real   (r8) :: flux_con_tol = 0.0_r8
  integer(in) :: flux_con_max_iter = 2

  !--- cold air outbreak parameters  (Mahrt & Sun 1995,MWR) -------------
  logical :: use_coldair_outbreak_mod = .false.
  real   (r8), parameter :: alpha  = 1.4_r8
  real   (r8), parameter :: maxscl = 2.0_r8     ! maximum wind scaling for flux
  real   (r8), parameter :: td0    = -10.0_r8   ! start t-ts for scaling

  character(*), parameter :: sourcefile = __FILE__

contains

  subroutine shr_flux_adjust_constants( &
       zvir, cpair, cpvir, karman, gravit, &
       latvap, latice, stebol, flux_convergence_tolerance, &
       flux_convergence_max_iteration, &
       coldair_outbreak_mod)

    ! Adjust local constants.  Used to support simple models.

    real(r8), optional, intent(in) :: zvir
    real(r8), optional, intent(in) :: cpair
    real(r8), optional, intent(in) :: cpvir
    real(r8), optional, intent(in) :: karman
    real(r8), optional, intent(in) :: gravit
    real(r8), optional, intent(in) :: latvap
    real(r8), optional, intent(in) :: latice
    real(r8), optional, intent(in) :: stebol
    real(r8), optional, intent(in)  :: flux_convergence_tolerance
    integer(in), optional, intent(in) :: flux_convergence_max_iteration
    logical, optional, intent(in) :: coldair_outbreak_mod
    !----------------------------------------------------------------------------

    if (present(zvir                          )) loc_zvir                 = zvir
    if (present(cpair                         )) loc_cpdair               = cpair
    if (present(cpvir                         )) loc_cpvir                = cpvir
    if (present(karman                        )) loc_karman               = karman
    if (present(gravit                        )) loc_g                    = gravit
    if (present(latvap                        )) loc_latvap               = latvap
    if (present(latice                        )) loc_latice               = latice
    if (present(stebol                        )) loc_stebol               = stebol
    if (present(flux_convergence_tolerance    )) flux_con_tol             = flux_convergence_tolerance
    if (present(flux_convergence_max_iteration)) flux_con_max_iter        = flux_convergence_max_iteration
    if (present(coldair_outbreak_mod          )) use_coldair_outbreak_mod = coldair_outbreak_mod

  end subroutine shr_flux_adjust_constants

  ! ----------------------------------------------------------------------------
  ! Description:
  !
  !     Internal atm/ocn flux calculation
  !
  ! Revision history:
  !     2002-Jun-10 - B. Kauffman - code migrated from cpl5 to cpl6
  !     2003-Apr-02 - B. Kauffman - taux & tauy now utilize ocn velocity
  !     2003-Apr-02 - B. Kauffman - tref,qref,duu10n mods as per Bill Large
  !     2006-Nov-07 - B. Kauffman - code migrated from cpl6 to share
  !
  !     2011-Mar-13 - J. Nusbaumer - Water Isotope ocean flux added.

  !     2019-May-16 - Jack Reeves Eyre (UA) and Kai Zhang (PNNL) -
  !                   Added COARE/Fairall surface flux scheme option
  !                   (ocn_surface_flux_scheme .eq. 1) based on code from
  !                   Thomas Toniazzo (Bjerknes Centre, Bergen) ”
  ! ----------------------------------------------------------------------------

  subroutine flux_atmOcn(logunit, nMax    , zbot    , ubot  , vbot  ,thbot   , &
                         qbot   , rainc   , s16O    , sHDO  , s18O  ,rbot    , &
                         tbot   , us      , vs      , pslv                   , &
                         ts     , mask    , seq_flux_atmocn_minwind          , &
                         sen    , lat     , lwup    , r16O  , rhdo  , r18O   , &
                         evap   , evap_16O, evap_HDO, evap_18O               , &
                         taux   , tauy    , tref    , qref                   , &
                         ocn_surface_flux_scheme    , add_gusts, duu10n      , &
                         ugust_out, u10res , ustar_sv, re_sv, ssq_sv, missval)

    use water_isotopes, only: wiso_flxoce

    integer    , intent(in ) :: logunit
    integer(in), intent(in ) :: nMax                     ! data vector length
    integer(in), intent(in ) :: mask     (nMax)          ! ocn domain mask       0 <=> out of domain
    integer(in), intent(in ) :: ocn_surface_flux_scheme
    logical    , intent(in ) :: add_gusts
    real   (r8), intent(in ) :: zbot     (nMax)          ! atm level height      (m)
    real   (r8), intent(in ) :: ubot     (nMax)          ! atm u wind            (m/s)
    real   (r8), intent(in ) :: vbot     (nMax)          ! atm v wind            (m/s)
    real   (r8), intent(in ) :: thbot    (nMax)          ! atm potential T       (K)
    real   (r8), intent(in ) :: qbot     (nMax)          ! atm specific humidity (kg/kg)
    real   (r8), intent(in ) :: rainc    (nMax)          ! atm precip for convective gustiness (kg/m^3) - RBN 24Nov2008/MDF 31Jan2022
    real   (r8), intent(in ) :: s16O     (nMax)          ! atm H216O tracer conc. (kg/kg)
    real   (r8), intent(in ) :: sHDO     (nMax)          ! atm HDO tracer conc.  (kg/kg)
    real   (r8), intent(in ) :: s18O     (nMax)          ! atm H218O tracer conc. (kg/kg)
    real   (r8), intent(in ) :: r16O     (nMax)          ! ocn H216O tracer ratio/Rstd
    real   (r8), intent(in ) :: rHDO     (nMax)          ! ocn HDO tracer ratio/Rstd
    real   (r8), intent(in ) :: r18O     (nMax)          ! ocn H218O tracer ratio/Rstd
    real   (r8), intent(in ) :: rbot     (nMax)          ! atm air density       (kg/m^3)
    real   (r8), intent(in ) :: tbot     (nMax)          ! atm T                 (K)
    real   (r8), intent(in ) :: pslv     (nMax)          ! atm sea level pressure(Pa)
    real   (r8), intent(in ) :: us       (nMax)          ! ocn u-velocity        (m/s)
    real   (r8), intent(in ) :: vs       (nMax)          ! ocn v-velocity        (m/s)
    real   (r8), intent(in ) :: ts       (nMax)          ! ocn temperature       (K)
    real   (r8), intent(in ) :: seq_flux_atmocn_minwind  ! minimum wind speed for atmocn (m/s)
    real   (r8), intent(out) :: sen      (nMax)          ! heat flux: sensible    (W/m^2)
    real   (r8), intent(out) :: lat      (nMax)          ! heat flux: latent      (W/m^2)
    real   (r8), intent(out) :: lwup     (nMax)          ! heat flux: lw upward   (W/m^2)
    real   (r8), intent(out) :: evap     (nMax)          ! water flux: evap  ((kg/s)/m^2)
    real   (r8), intent(out) :: evap_16O (nMax)          ! water flux: evap ((kg/s/m^2)
    real   (r8), intent(out) :: evap_HDO (nMax)          ! water flux: evap ((kg/s)/m^2)
    real   (r8), intent(out) :: evap_18O (nMax)          ! water flux: evap ((kg/s/m^2)
    real   (r8), intent(out) :: taux     (nMax)          ! surface stress, zonal      (N)
    real   (r8), intent(out) :: tauy     (nMax)          ! surface stress, maridional (N)
    real   (r8), intent(out) :: tref     (nMax)          ! diag:  2m ref height T     (K)
    real   (r8), intent(out) :: qref     (nMax)          ! diag:  2m ref humidity (kg/kg)
    real   (r8), intent(out) :: duu10n   (nMax)          ! diag: 10m wind speed squared (m/s)^2
    real   (r8), intent(out) :: ugust_out(nMax)          ! diag: gustiness addition to U10 (m/s)
    real   (r8), intent(out) :: u10res   (nMax)          ! diag: gustiness addition to U10 (m/s)
    real   (r8), intent(out), optional :: ustar_sv(nMax) ! diag: ustar
    real   (r8), intent(out), optional :: re_sv   (nMax) ! diag: sqrt of exchange coefficient (water)
    real   (r8), intent(out), optional :: ssq_sv  (nMax) ! diag: sea surface humidity  (kg/kg)
    real   (r8), intent(in ), optional :: missval        ! masked value

    real   (r8), parameter :: zref  = 10.0_r8 ! reference height           (m)
    real   (r8), parameter :: ztref =  2.0_r8 ! reference height for air T (m)
    !!++ Large only
    !real(r8),parameter :: cexcd  = 0.0346_r8 ! ratio Ch(water)/CD
    !real(r8),parameter :: chxcds = 0.018_r8  ! ratio Ch(heat)/CD for stable case
    !real(r8),parameter :: chxcdu = 0.0327_r8 ! ratio Ch(heat)/CD for unstable case
    !!++ COARE only
    real   (r8), parameter :: zpbl = 700.0_r8 ! PBL depth [m] for gustiness parametriz.

    integer(in) n, iter
    real(r8)    :: vmag               ! surface wind magnitude   (m/s)
    real(r8)    :: ssq                ! sea surface humidity     (kg/kg)
    real(r8)    :: delt               ! potential T difference   (K)
    real(r8)    :: delq               ! humidity difference      (kg/kg)
    real(r8)    :: stable             ! stability factor
    real(r8)    :: rdn                ! sqrt of neutral exchange coeff (momentum)
    real(r8)    :: rhn                ! sqrt of neutral exchange coeff (heat)
    real(r8)    :: ren                ! sqrt of neutral exchange coeff (water)
    real(r8)    :: rd                 ! sqrt of exchange coefficient (momentum)
    real(r8)    :: rh                 ! sqrt of exchange coefficient (heat)
    real(r8)    :: re                 ! sqrt of exchange coefficient (water)
    real(r8)    :: ustar              ! ustar
    real(r8)    :: ustar_prev
    real(r8)    :: qstar              ! qstar
    real(r8)    :: tstar              ! tstar
    real(r8)    :: hol                ! H (at zbot) over L
    real(r8)    :: xsq                ! ?
    real(r8)    :: xqq                ! ?
    !!++ Large only
    real(r8)    :: psimh              ! stability function at zbot (momentum)
    real(r8)    :: psixh              ! stability function at zbot (heat and water)
    real(r8)    :: psix2              ! stability function at ztref reference height
    real(r8)    :: alz                ! ln(zbot/zref)
    real(r8)    :: al2                ! ln(zref/ztref)
    real(r8)    :: u10n               ! 10m neutral wind
    real(r8)    :: tau                ! stress at zbot
    real(r8)    :: cp                 ! specific heat of moist air
    real(r8)    :: fac                ! vertical interpolation factor
    real(r8)    :: spval              ! local missing value
    real(r8)    :: wind0              ! resolved large-scale 10m wind (no gust added)
    !!++ COARE only
    real(r8)    :: zo, zot, zoq       ! roughness lengths
    real(r8)    :: hsb, hlb           ! sens & lat heat flxs at zbot
    real(r8)    :: trf, qrf, urf, vrf ! reference-height quantities


    !--- local functions --------------------------------
    real(r8)    :: qsat               ! function: the saturation humididty of air (kg/m^3)
    !!++ Large only (formula v*=[c4/U10+c5+c6*U10]*U10 in Large et al. 1994)
    real(r8)    :: cdn                ! function: neutral drag coeff at 10m
    !!++ Large only (stability functions)
    real(r8)    :: psimhu             ! function: unstable part of psimh
    real(r8)    :: psixhu             ! function: unstable part of psimx
    real(r8)    :: Umps               ! dummy arg ~ wind velocity (m/s)
    real(r8)    :: Tk                 ! dummy arg ~ temperature (K)
    real(r8)    :: xd                 ! dummy arg ~ ?
    !--- for cold air outbreak calc --------------------------------
    real(r8)    :: tdiff(nMax)        ! tbot - ts
    real(r8)    :: vscl
    real(r8)    :: ugust              ! function: gustiness as a function of convective rainfall.  
    real(r8)    :: gprec              ! convective rainfall argument for ugust

    qsat(Tk) = 640380.0_r8 / exp(5107.4_r8 / Tk)

    ! Large and Yeager 2009
    cdn(Umps) = 0.0027_r8 / min(33.0000_r8, Umps) + 0.000142_r8 + &
         0.0000764_r8 * min(33.0000_r8, Umps) - 3.14807e-13_r8 * min(33.0000_r8, Umps)**6
    ! Capped Large and Pond by wind
    ! cdn(Umps) = 0.0027_r8 / min(30.0_r8,Umps) + 0.000142_r8 + 0.0000764_r8 * min(30.0_r8,Umps)
    ! Capped Large and Pond by Cd
    ! cdn(Umps) = min(0.0025_r8, (0.0027_r8 / Umps + 0.000142_r8 + 0.0000764_r8 * Umps ))
    ! Large and Pond
    ! cdn(Umps) = 0.0027_r8 / Umps + 0.000142_r8 + 0.0000764_r8 * Umps

    psimhu(xd) = log((1.0_r8 + xd * (2.0_r8 + xd)) * (1.0_r8 + xd * xd) / 8.0_r8) - 2.0_r8 * atan(xd) + 1.571_r8
    psixhu(xd) = 2.0_r8 * log((1.0_r8 + xd * xd)/2.0_r8)

    ! Convective gustiness appropriate for input precipitation.
    ! Following Regelsperger et al. (2000, J. Clim)
    ! Ug = log(1.0+6.69R-0.476R^2)
    ! Coefficients X by 8640 for mm/s (from cam) -> cm/day (for above forumla)
    ugust(gprec) = log(1.0_r8 + 57801.6_r8 * gprec - 3.55332096e7_r8 * (gprec**2))

    character(*), parameter :: subName = '(flux_atmOcn) '
    character(*), parameter ::     F00 = "('(flux_atmOcn) ',4a)"

    !-------------------------------------------------------------------------------
    ! PURPOSE:
    !   computes atm/ocn surface fluxes
    !
    ! NOTES:
    !   o all fluxes are positive downward
    !   o net heat flux = net sw + lw up + lw down + sen + lat
    !   o here, tstar = <WT>/U*, and qstar = <WQ>/U*.
    !   o wind speeds should all be above a minimum speed (eg. 1.0 m/s)
    !
    ! ASSUMPTIONS:
    !  Large:
    !   o Neutral 10m drag coeff: cdn = .0027/U10 + .000142 + .0000764 U10
    !   o Neutral 10m stanton number: ctn = .0327 sqrt(cdn), unstable
    !                                 ctn = .0180 sqrt(cdn), stable
    !   o Neutral 10m dalton number:  cen = .0346 sqrt(cdn)
    !   o The saturation humidity of air at T(K): qsat(T)  (kg/m^3)
    !  COARE:
    !   o use COAREv3.0 function (tht 22/11/2013)
    !-------------------------------------------------------------------------------

    if (debug > 0) write(logunit, F00) 'enter'

    if (present(missval)) then
      spval = missval
    else
      spval = shr_const_spval
    end if
    u10n  = spval
    rh    = spval
    psixh = spval
    hol   = spval

    !--- for cold air outbreak calc --------------------------------
    tdiff = tbot - ts

    !.................................................................
    ! ocn_surface_flux_scheme = 0 : Default CESM1.2
    !                         = 1 : COARE algorithm
    !                         = 2 : UA algorithm (separate subroutine)
    !.................................................................

    ! Default flux scheme.
    if (ocn_surface_flux_scheme == 0) then
      al2 = log(zref / ztref)
      do n = 1, nMax
        if (mask(n) /= 0) then
          !--- compute some needed quantities ---
          if (add_gusts) then 
            vmag = max(seq_flux_atmocn_minwind, sqrt((ubot(n) - us(n))**2 + (vbot(n) - vs(n))**2 + (1.0_r8 * ugust(min(rainc(n), 6.94444e-4_r8))**2)))
            ugust_out(n) = ugust(min(rainc(n), 6.94444e-4_r8))
          else 
            vmag = max(seq_flux_atmocn_minwind, sqrt((ubot(n) - us(n))**2 + (vbot(n) - vs(n))**2))
            ugust_out(n) = 0
          end if
          wind0 = max(seq_flux_atmocn_minwind, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) )
          if (use_coldair_outbreak_mod) then
            ! Cold Air Outbreak Modification:
            ! Increase windspeed for negative tbot-ts
            ! based on Mahrt & Sun 1995, MWR
            if (tdiff(n) < td0) then
              ! if add_gusts wind0 and vmag are different, both need this factor.
              vscl = min((1.0_r8 + alpha * (abs(tdiff(n) - td0)**0.5_r8 / abs(vmag))), maxscl)
              vmag = vmag * vscl
              vscl = min((1.0_r8 + alpha * (abs(tdiff(n) - td0)**0.5_r8 / abs(wind0))),maxscl)
              wind0 = wind0 * vscl
            end if
          end if
          ssq    = 0.98_r8 * qsat(ts(n)) / rbot(n)   ! sea surf hum (kg/kg)
          delt   = thbot(n) - ts(n)                  ! pot temp diff (K)
          delq   = qbot(n) - ssq                     ! spec hum dif (kg/kg)
          alz    = log(zbot(n) / zref)
          cp     = loc_cpdair * (1.0_r8 + loc_cpvir * ssq)

          !------------------------------------------------------------
          ! first estimate of Z/L and ustar, tstar and qstar
          !------------------------------------------------------------
          !--- neutral coefficients, z/L = 0.0 ---
          stable = 0.5_r8 + sign(0.5_r8 , delt)
          rdn    = sqrt(cdn(vmag))
          rhn    = (1.0_r8-stable) * 0.0327_r8 + stable * 0.018_r8
                 ! (1.0_r8-stable) * chxcdu + stable * chxcds
          ren    = 0.0346_r8 ! cexcd

          !--- ustar, tstar, qstar ---
          ustar = rdn * vmag
          tstar = rhn * delt
          qstar = ren * delq
          ustar_prev = ustar*2.0_r8
          iter = 0
          do while (abs((ustar - ustar_prev)/ustar) > flux_con_tol .and. iter < flux_con_max_iter)
            iter = iter + 1
            ustar_prev = ustar
            !--- compute stability & evaluate all stability functions ---
            hol    = loc_karman * loc_g * zbot(n) * (tstar / thbot(n) + qstar / (1.0_r8 / loc_zvir+qbot(n))) / ustar**2
            hol    = sign(min(abs(hol), 10.0_r8), hol)
            stable = 0.5_r8 + sign(0.5_r8 , hol)
            xsq    = max(sqrt(abs(1.0_r8 - 16.0_r8*hol)) , 1.0_r8)
            xqq    = sqrt(xsq)
            psimh  = -5.0_r8 * hol * stable + (1.0_r8 - stable) * psimhu(xqq)
            psixh  = -5.0_r8 * hol * stable + (1.0_r8 - stable) * psixhu(xqq)

            !--- shift wind speed using old coefficient ---
            rd     = rdn / (1.0_r8 + max(rdn / loc_karman * (alz - psimh), -0.5_r8))
            u10n   = vmag * rd / rdn

            !--- update transfer coeffs at 10m and neutral stability ---
            rdn    = sqrt(cdn(u10n))
            ren    = 0.0346_r8 ! cexcd
            rhn    = (1.0_r8 - stable) * 0.0327_r8 + stable * 0.018_r8
                   ! (1.0_r8 - stable) * chxcdu + stable * chxcds

            !--- shift all coeffs to measurement height and stability ---
            rd     = rdn / (1.0_r8 + rdn / loc_karman * (alz - psimh))
            rh     = rhn / (1.0_r8 + rhn / loc_karman * (alz - psixh))
            re     = ren / (1.0_r8 + ren / loc_karman * (alz - psixh))

            !--- update ustar, tstar, qstar using updated, shifted coeffs --
            ustar  = rd * vmag
            tstar  = rh * delt
            qstar  = re * delq
          end do
          if (iter < 1) then
            write(logunit, *) ustar,ustar_prev,flux_con_tol,flux_con_max_iter
            call shr_sys_abort('No iterations performed in flux_atmocn_mod')
          end if
          !------------------------------------------------------------
          ! compute the fluxes
          !------------------------------------------------------------

          tau = rbot(n) * ustar * ustar

          !--- momentum flux ---
          taux(n) = tau * (ubot(n) - us(n)) / vmag
          tauy(n) = tau * (vbot(n) - vs(n)) / vmag

          !--- heat flux ---
          sen (n) =          cp * tau * tstar / ustar
          lat (n) =  loc_latvap * tau * qstar / ustar
          lwup(n) = -loc_stebol * ts(n)**4

          !--- water flux ---
          evap(n) = lat(n)/loc_latvap

          !---water isotope flux ---
          call wiso_flxoce(2, rbot(n), zbot(n), s16O(n), ts(n), r16O(n), ustar, re, ssq, evap_16O(n), qbot(n), evap(n))
          call wiso_flxoce(3, rbot(n), zbot(n), sHDO(n), ts(n), rHDO(n), ustar, re, ssq, evap_HDO(n), qbot(n), evap(n))
          call wiso_flxoce(4, rbot(n), zbot(n), s18O(n), ts(n), r18O(n), ustar, re, ssq, evap_18O(n), qbot(n), evap(n))

          !------------------------------------------------------------
          ! compute diagnositcs: 2m ref T & Q, 10m wind speed squared
          !------------------------------------------------------------
          hol     = hol * ztref / zbot(n)
          xsq     = max(1.0_r8, sqrt(abs(1.0_r8 - 16.0_r8 * hol)))
          xqq     = sqrt(xsq)
          psix2   = -5.0_r8*hol*stable + (1.0_r8-stable)*psixhu(xqq)
          fac     = (rh/loc_karman) * (alz + al2 - psixh + psix2 )
          tref(n) = thbot(n) - delt*fac
          tref(n) = tref(n) - 0.01_r8 * ztref   ! pot temp to temp correction
          fac     = (re / loc_karman) * (alz + al2 - psixh + psix2 )
          qref(n) =  qbot(n) - delq * fac

          duu10n(n) = u10n * u10n            ! 10m wind speed squared
          u10res(n) = u10n * (wind0 / vmag)  ! resolved 10m wind

          !------------------------------------------------------------
          ! optional diagnostics, needed for water tracer fluxes (dcn)
          !------------------------------------------------------------
          if (present(ustar_sv)) ustar_sv(n) = ustar
          if (present(re_sv   )) re_sv(n)    = re
          if (present(ssq_sv  )) ssq_sv(n)   = ssq
        else
          !------------------------------------------------------------
          ! no valid data here -- out of domain
          !------------------------------------------------------------
          sen      (n) = spval  ! sensible         heat flux  (W/m^2)
          lat      (n) = spval  ! latent           heat flux  (W/m^2)
          lwup     (n) = spval  ! long-wave upward heat flux  (W/m^2)
          evap     (n) = spval  ! evaporative water flux ((kg/s)/m^2)
          evap_16O (n) = spval  ! water tracer flux (kg/s)/m^2)
          evap_HDO (n) = spval  ! HDO tracer flux  (kg/s)/m^2)
          evap_18O (n) = spval  ! H218O tracer flux (kg/s)/m^2)
          taux     (n) = spval  ! x surface stress (N)
          tauy     (n) = spval  ! y surface stress (N)
          tref     (n) = spval  ! 2m reference height temperature (K)
          qref     (n) = spval  ! 2m reference height humidity (kg/kg)
          duu10n   (n) = spval  ! 10m wind speed squared (m/s)^2
          ugust_out(n) = spval  ! gustiness addition (m/s)
          u10res   (n) = spval  ! 10m resolved wind (no gusts) (m/s) 

          if (present(ustar_sv)) ustar_sv(n) = spval
          if (present(re_sv   )) re_sv   (n) = spval
          if (present(ssq_sv  )) ssq_sv  (n) = spval
        end if
      end do
    else if (ocn_surface_flux_scheme == 1) then
      !.................................
      ! Use COARE algorithm
      !.................................
      do n = 1, nMax
        if (mask(n) /= 0) then
          !--- compute some needed quantities ---
          vmag = max(seq_flux_atmocn_minwind, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) )

          if (use_coldair_outbreak_mod) then
            ! Cold Air Outbreak Modification:
            ! Increase windspeed for negative tbot-ts
            ! based on Mahrt & Sun 1995,MWR
            if (tdiff(n) < td0) then
              vscl = min((1 + alpha * (abs(tdiff(n) - td0)**0.5_r8 / abs(vmag))), maxscl)
              vmag = vmag * vscl
            end if
          end if
          ssq = 0.98_r8 * qsat(ts(n)) / rbot(n)   ! sea surf hum (kg/kg)

          call cor30a(ubot(n), vbot(n), tbot(n), qbot(n), rbot(n), &               ! in atm params
                      us  (n), vs  (n), ts  (n), ssq             , &               ! in surf params
                      zpbl   , zbot(n), zbot(n), zref   , ztref  , ztref, &        ! in heights
                      tau    , hsb    , hlb                      , &               ! out: fluxes
                      zo     , zot    , zoq    , hol    , ustar  , tstar, qstar, & ! out: ss scales
                      rd     , rh     , re                       , &               ! out: exch. coeffs
                      trf    , qrf    , urf    , vrf)                              ! out: reference-height params

          ! For the sake of maintaining same defs
          hol = zbot(n) / hol
          rd  = sqrt(rd)
          rh  = sqrt(rh)
          re  = sqrt(re)

          !--- momentum flux ---
          taux(n) = tau * (ubot(n) - us(n)) / vmag
          tauy(n) = tau * (vbot(n) - vs(n)) / vmag

          !--- heat flux ---
          sen (n) =  hsb
          lat (n) =  hlb
          lwup(n) = -shr_const_stebol * ts(n)**4

          !--- water flux ---
          evap(n) = lat(n) / shr_const_latvap

          !---water isotope flux ---
          call wiso_flxoce(2, rbot(n), zbot(n), s16O(n), ts(n), r16O(n), ustar, re, ssq, evap_16O(n), qbot(n), evap(n))
          call wiso_flxoce(3, rbot(n), zbot(n), sHDO(n), ts(n), rHDO(n), ustar, re, ssq, evap_HDO(n), qbot(n), evap(n))
          call wiso_flxoce(4, rbot(n), zbot(n), s18O(n), ts(n), r18O(n), ustar, re, ssq, evap_18O(n), qbot(n), evap(n))

          !------------------------------------------------------------
          ! compute diagnositcs: 2m ref T & Q, 10m wind speed squared
          !------------------------------------------------------------
          tref(n) = trf
          qref(n) = qrf
          duu10n(n) = urf**2 + vrf**2

          !------------------------------------------------------------
          ! optional diagnostics, needed for water tracer fluxes (dcn)
          !------------------------------------------------------------
          if (present(ustar_sv)) ustar_sv(n) = ustar
          if (present(re_sv   )) re_sv(n)    = re
          if (present(ssq_sv  )) ssq_sv(n)   = ssq
        else
          !------------------------------------------------------------
          ! no valid data here -- out of domain
          !------------------------------------------------------------
          sen      (n) = spval  ! sensible         heat flux  (W/m^2)
          lat      (n) = spval  ! latent           heat flux  (W/m^2)
          lwup     (n) = spval  ! long-wave upward heat flux  (W/m^2)
          evap     (n) = spval  ! evaporative water flux ((kg/s)/m^2)
          evap_16O (n) = spval  ! water tracer flux (kg/s)/m^2)
          evap_HDO (n) = spval  ! HDO tracer flux  (kg/s)/m^2)
          evap_18O (n) = spval  ! H218O tracer flux (kg/s)/m^2)
          taux     (n) = spval  ! x surface stress (N)
          tauy     (n) = spval  ! y surface stress (N)
          tref     (n) = spval  !  2m reference height temperature (K)
          qref     (n) = spval  !  2m reference height humidity (kg/kg)
          duu10n   (n) = spval  ! 10m wind speed squared (m/s)^2

          if (present(ustar_sv)) ustar_sv(n) = spval
          if (present(re_sv   )) re_sv   (n) = spval
          if (present(ssq_sv  )) ssq_sv  (n) = spval
        end if
      end do
    else if (ocn_surface_flux_scheme == 2) then
      call flux_atmOcn_UA(logunit, &
                          nMax, zbot, ubot, vbot, thbot, &
                          qbot, s16O, sHDO, s18O, rbot , &
                          tbot, pslv, us  , vs  , ts   , &
                          mask, sen , lat , lwup, r16O , &
                          rhdo, r18O, evap, evap_16O, evap_HDO, evap_18O, &
                          taux, tauy, tref, qref, duu10n, ustar_sv, re_sv, &
                          ssq_sv, missval)
    else
      call shr_sys_abort(subName // ' subroutine flux_atmOcn requires ocn_surface_flux_scheme = 0, 1 or 2')
    end if

  end subroutine flux_atmocn

  !===============================================================================
  ! !IROUTINE: flux_atmOcn_UA -- internal atm/ocn flux calculation
  !
  ! !DESCRIPTION:
  !
  !     Internal atm/ocn flux calculation
  !     using University of Arizona method.
  !
  !     Reference:
  !         Zeng, X., M. Zhao, and R.E. Dickinson, 1998: Intercomparison of Bulk
  !             Aerodynamic Algorithms for the Computation of Sea Surface Fluxes
  !             Using TOGA COARE and TAO Data. J. Climate, 11, 2628–2644,
  !             https://doi.org/10.1175/1520-0442(1998)011<2628%3AIOBAAF>2.0.CO%3B2
  !
  !     Equation numbers are from this paper.
  !
  ! !REVISION HISTORY:
  !     2017-Aug-28 - J. Reeves Eyre - code re-written for E3SM
  !     2018-Oct-30 - J. Reeves Eyre - bug fix and add
  !                   convective gustiness.
  !     2019-May-08 - J. Reeves Eyre - remove convective gustiness
  !                   and add cold air outbreak modification.
  !===============================================================================
  SUBROUTINE flux_atmOcn_UA(logunit,  &
       &               nMax  ,zbot  ,ubot  ,vbot  ,thbot ,  &
       &               qbot  ,s16O  ,sHDO  ,s18O  ,rbot  ,   &
       &               tbot  , pslv ,us    , vs   ,   &
       &               ts    ,mask  ,sen   ,lat   ,lwup  ,   &
       &               r16O, rhdo, r18O, &
       &               evap  ,evap_16O, evap_HDO, evap_18O, &
       &               taux  ,tauy  ,tref  ,qref  ,   &
       &               duu10n,  ustar_sv   ,re_sv ,ssq_sv,   &
       &               missval)


    ! !USES:
    use water_isotopes, only: wiso_flxoce !subroutine used to calculate water isotope fluxes.

    implicit none

    ! !INPUT/OUTPUT PARAMETERS:

    !--- input arguments --------------------------------
    integer    ,intent(in) :: logunit
    integer    ,intent(in) :: nMax        ! data vector length
    integer    ,intent(in) :: mask (nMax) ! ocn domain mask       0 <=> out of domain
    real(r8)   ,intent(in) :: zbot (nMax) ! atm level height      (m)
    real(r8)   ,intent(in) :: ubot (nMax) ! atm u wind            (m/s)
    real(r8)   ,intent(in) :: vbot (nMax) ! atm v wind            (m/s)
    real(r8)   ,intent(in) :: thbot(nMax) ! atm potential T       (K)
    real(r8)   ,intent(in) :: qbot (nMax) ! atm specific humidity (kg/kg)
    real(r8)   ,intent(in) :: s16O (nMax) ! atm H216O tracer conc. (kg/kg)
    real(r8)   ,intent(in) :: sHDO (nMax) ! atm HDO tracer conc.  (kg/kg)
    real(r8)   ,intent(in) :: s18O (nMax) ! atm H218O tracer conc. (kg/kg)
    real(r8)   ,intent(in) :: r16O (nMax) ! ocn H216O tracer ratio/Rstd
    real(r8)   ,intent(in) :: rHDO (nMax) ! ocn HDO tracer ratio/Rstd
    real(r8)   ,intent(in) :: r18O (nMax) ! ocn H218O tracer ratio/Rstd
    real(r8)   ,intent(in) :: rbot (nMax) ! atm air density       (kg/m^3)
    real(r8)   ,intent(in) :: tbot (nMax) ! atm T                 (K)
    real(r8)   ,intent(in) :: pslv (nMax) ! sea level pressure    (Pa)
    real(r8)   ,intent(in) :: us   (nMax) ! ocn u-velocity        (m/s)
    real(r8)   ,intent(in) :: vs   (nMax) ! ocn v-velocity        (m/s)
    real(r8)   ,intent(in) :: ts   (nMax) ! ocn temperature       (K)

    !--- output arguments -------------------------------
    real(r8),intent(out)  ::  sen  (nMax) ! heat flux: sensible    (W/m^2)
    real(r8),intent(out)  ::  lat  (nMax) ! heat flux: latent      (W/m^2)
    real(r8),intent(out)  ::  lwup (nMax) ! heat flux: lw upward   (W/m^2)
    real(r8),intent(out)  ::  evap (nMax) ! water flux: evap  ((kg/s)/m^2)
    real(r8),intent(out)  ::  evap_16O (nMax) ! water flux: evap ((kg/s/m^2)
    real(r8),intent(out)  ::  evap_HDO (nMax) ! water flux: evap ((kg/s)/m^2)
    real(r8),intent(out)  ::  evap_18O (nMax) ! water flux: evap ((kg/s/m^2)
    real(r8),intent(out)  ::  taux (nMax) ! surface stress, zonal      (N)
    real(r8),intent(out)  ::  tauy (nMax) ! surface stress, maridional (N)
    real(r8),intent(out)  ::  tref (nMax) ! diag:  2m ref height T     (K)
    real(r8),intent(out)  ::  qref (nMax) ! diag:  2m ref humidity (kg/kg)
    real(r8),intent(out)  :: duu10n(nMax) ! diag: 10m wind speed squared (m/s)^2

    real(r8),intent(out),optional :: ustar_sv(nMax) ! diag: ustar
    real(r8),intent(out),optional :: re_sv   (nMax) ! diag: sqrt of exchange coefficient (water)
    real(r8),intent(out),optional :: ssq_sv  (nMax) ! diag: sea surface humidity  (kg/kg)

    real(r8),intent(in) ,optional :: missval        ! masked value

    !--- local constants --------------------------------
    real(r8),parameter :: zetam = -1.574_r8 ! Very unstable zeta cutoff for momentum (-)
    real(r8),parameter :: zetat = -0.465_r8 ! Very unstable zeta cutoff for T/q (-)
    real(r8),parameter :: umin  = 0.1_r8    ! minimum wind speed       (m/s)
    real(r8),parameter :: zref  = 10.0_r8   ! reference height           (m)
    real(r8),parameter :: ztref = 2.0_r8    ! reference height for air T (m)
    real(r8),parameter :: beta = 1.0_r8     ! constant used in W* calculation (-)
    real(r8),parameter :: zpbl = 1000.0_r8  ! PBL height used in W* calculation (m)
    real(r8),parameter :: gamma = 0.0098_r8         ! Dry adiabatic lapse rate (K/m)
    real(r8),parameter :: onethird = 1.0_r8/3.0_r8  ! Used repeatedly.

    !--- local variables --------------------------------
    integer(in) :: n          ! vector loop index
    integer(in) :: i          ! iteration loop index
    real(r8)    :: vmag_abs   ! surface wind magnitude   (m s-1)
    real(r8)    :: vmag_rel   ! surface wind magnitude relative to
    ! surface current   (m s-1)
    real(r8)    :: vmag       ! surface wind magnitude with large
    ! eddy correction and minimum value (m s-1)
    ! (This can change on each iteration.)
    real(r8)    :: thv        ! virtual temperature      (K)
    real(r8)    :: ssq        ! sea surface humidity     (kg/kg)
    real(r8)    :: delth      ! potential T difference   (K)
    real(r8)    :: delthv     ! virtual potential T difference   (K)
    real(r8)    :: delq       ! humidity difference      (kg/kg)
    real(r8)    :: ustar      ! friction velocity (m s-1)
    real(r8)    :: qstar      ! humidity scaling parameter (kg/kg)
    real(r8)    :: tstar      ! temperature scaling parameter (K)
    real(r8)    :: thvstar    ! virtual temperature scaling parameter (K)
    real(r8)    :: wstar      ! convective velocity scale (m s-1)
    real(r8)    :: zeta       ! dimensionless height (z / Obukhov length)
    real(r8)    :: obu        ! Obukhov length (m)
    real(r8)    :: tau        ! magnitude of wind stress (N m-2)
    real(r8)    :: cp         ! specific heat of moist air (J kg-1 K-1)
    real(r8)    :: xlv        ! Latent heat of vaporization (J kg-1)
    real(r8)    :: visa       ! Kinematic viscosity of dry air (m2 s-1)
    real(r8)    :: tbot_oC    ! Temperature used in visa (deg C)
    real(r8)    :: rb         ! Bulk Richardson number (-)
    real(r8)    :: zo         ! Roughness length for momentum (m)
    real(r8)    :: zoq        ! Roughness length for moisture (m)
    real(r8)    :: zot        ! Roughness length for heat (m)
    real(r8)    :: u10        ! 10-metre wind speed (m s-1)
    real(r8)    :: re         ! Moisture exchange coefficient for compatibility
    ! with default algorithm.
    real(r8)    :: spval      ! local missing value
    real(r8)    :: loc_epsilon  ! Ratio of gas constants (-)

    !--- for cold air outbreak calc --------------------------------
    real(r8)    :: tdiff(nMax)  ! tbot - ts
    real(r8)    :: vscl

    !--- formats ----------------------------------------
    character(*),parameter :: subName = '(flux_atmOcn) '
    character(*),parameter ::   F00 = "('(flux_atmOcn) ',4a)"

    !-----
    ! Straight from original subroutine.
    if (debug > 0) write(logunit,F00) "enter"

    if (present(missval)) then
       spval = missval
    else
       spval = shr_const_spval
    endif
    !-----

    ! Evaluate loc_epsilon.
    loc_epsilon = 1.0_r8 / (1.0_r8 + loc_zvir)

    !--- for cold air outbreak calc --------------------------------
    tdiff = tbot - ts

    ! Loop over grid points.
    DO n=1,nMax
       if (mask(n) /= 0) then

          !-----Calculate some required near surface variables.---------
          vmag_abs = sqrt( ubot(n)**2 + vbot(n)**2 )
          vmag_rel = sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2 )

          ! For Cold Air Outbreak Modification (based on Mahrt & Sun 1995,MWR):
          if (use_coldair_outbreak_mod) then
             ! Increase windspeed for negative tbot-ts
             if (tdiff(n).lt.td0) then
                vscl=min((1._r8+alpha*(abs(tdiff(n)-td0)**0.5_r8/abs(vmag_rel))),maxscl)
                vmag_rel=vmag_rel*vscl
             endif
          endif

          delth = thbot(n) - ts(n)                     ! Pot. temp. difference with surface (K)
          ! Note this is equivalent to Zeng et al
          ! (1998) version = delt + 0.0098*zbot
          thv = thbot(n)*(1.0_r8+0.61_r8*qbot(n))      ! Virtual potential temperature (K)
          ! EQN (17):
          !ssq = 0.98_r8 * qsat_ua(ts(n),ps, &          ! Surface specific humidity (kg kg-1)
          !                        loc_epsilon)
          ssq = 0.98_r8 * qsat_ua(ts(n),pslv(n), &     ! Surface specific humidity (kg kg-1)
               loc_epsilon)
          delq = qbot(n) - ssq                         ! Difference to surface (kg kg-1)
          delthv = delth*(1.0_r8+0.61_r8*qbot(n)) + &  ! Difference of virtual potential
               & 0.61_r8*thbot(n)*delq               ! temperature with surface (K)

          xlv = 1.0e+6_r8 * &                           ! Latent heat of vaporization (J kg-1)
               & (2.501_r8 - 0.00237_r8 * (ts(n) - loc_tkfrz))
          tbot_oC = tbot(n) - loc_tkfrz
          visa = 1.326e-5_r8 * (1.0_r8 + &             ! Kinematic viscosity of dry
               & 6.542e-3_r8*tbot_oC + &               ! air (m2 s-1) from Andreas (1989)
               & 8.301e-6_r8*tbot_oC*tbot_oC - &       ! CRREL Rep. 89-11
               & 4.84e-9_r8*tbot_oC*tbot_oC*tbot_oC)
          cp = loc_cpdair*(1.0_r8 + loc_cpvir*ssq)     ! specific heat of moist air (J kg-1 K-1)

          !-----Initial values of u* and convective velocity.-----------
          ustar = 0.06_r8
          wstar = 0.5_r8
          ! Update wind speed if unstable regime.
          if (delthv.lt.0.0_r8) then
             ! EQN (19)
             vmag = sqrt( vmag_rel**2 + beta*beta*wstar*wstar )
          else
             ! EQN (18)
             vmag = max(umin,vmag_rel)
          endif

          !-----Iterate to compute new u* and z0.-----------------------
          do i = 1,5
             ! EQN (24)
             zo = 0.013_r8*ustar*ustar/loc_g + 0.11_r8*visa/ustar
             ! EQN (9) assuming neutral
             ustar = loc_karman*vmag/log(zbot(n)/zo)
          enddo

          !-----Assess stability.---------------------------------------
          rb = loc_g*zbot(n)*delthv / (thv*vmag*vmag)    ! bulk Richardson number

          if(rb.ge.0.0_r8) then
             ! Neutral or stable: EQNs (4), (9), (13) and definition of rb.
             zeta = rb*log(zbot(n)/zo) / &
                  & (1.0_r8 - 5.0_r8*min(rb,0.19_r8))
          else
             ! Unstable: EQNs (4), (8), (12) and definition of rb.
             zeta = rb*log(zbot(n)/zo)
          endif

          obu = zbot(n)/zeta                             ! Obukhov length
          obu = sign(max(zbot(n)/10.0_r8, abs(obu)), obu)

          !-----Main iterations (2-10 iterations would be fine).-------
          do i=1,10

             ! Update roughness lengths.
             call rough_ua(zo,zot,zoq,ustar,visa)

             ! Wind variables.
             zeta = zbot(n) / obu
             if (zeta.lt.zetam) then
                ! Very unstable regime
                ! EQN (7) with extra z0 term.
                ustar = loc_karman * vmag / (log(zetam*obu/zo) - &
                     & psi_ua(1_IN, zetam) + &
                     & psi_ua(1_IN, zo/obu) + &
                     & 1.14_r8 * ((-zeta)**onethird - (-zetam)**onethird) )
             else if (zeta.lt.0.0_r8) then
                ! Unstable regime
                ! EQN (8) with extra z0 term.
                ustar = loc_karman * vmag / (log(zbot(n)/zo) - &
                     & psi_ua(1_IN,zeta) + psi_ua(1_IN,zo/obu) )
             else if (zeta.le.1.0_r8) then
                ! Stable regime
                ! EQN (9) with extra z0 term.
                ustar = loc_karman * vmag / (log(zbot(n)/zo) + &
                     & 5.0_r8*zeta - 5.0_r8*zo/obu)
             else
                ! Very stable regime
                ! EQN (10) with extra z0 term.
                ustar = loc_karman * vmag / (log(obu/zo) + 5.0_r8 - &
                     &  5.0_r8*zo/obu + &
                     &  (5.0_r8*log(zeta) + zeta - 1.0_r8) )
             endif

             ! Temperature variables.
             if(zeta.lt.zetat) then
                ! Very unstable regime
                ! EQN (11) with extra z0 term.
                tstar = loc_karman * delth / (log(zetat*obu/zot) - &
                     & psi_ua(2_IN, zetat) + &
                     & psi_ua(2_IN, zot/obu) + &
                     & 0.8_r8*((-zetat)**(-onethird) - (-zeta)**(-onethird)) )
             else if (zeta.lt.0.0_r8) then
                ! Unstable regime
                ! EQN (12) with extra z0 term.
                tstar = loc_karman * delth / &
                     & (log(zbot(n)/zot) - psi_ua(2_IN,zeta) + psi_ua(2_IN,zot/obu))
             else if (zeta.le.1.0_r8) then
                ! Stable regime
                ! EQN (13) with extra z0 term.
                tstar = loc_karman * delth / (log(zbot(n)/zot) + &
                     &   5.0_r8*zeta - 5.0_r8*zot/obu)
             else
                ! Very stable regime
                ! EQN (14) with extra z0 term.
                tstar = loc_karman * delth / (log(obu/zot) + &
                     &   5.0_r8 - 5.0_r8*zot/obu  + &
                     &   (5.0_r8*log(zeta) + zeta - 1.0_r8) )
             endif

             ! Humidity variables.
             ! This is done with re to give variable to save out like
             ! in old algorithm.
             if (zeta.lt.zetat) then
                ! Very unstable regime
                ! EQN (11) with extra z0 term.
                re = loc_karman / (log(zetat*obu/zoq) - psi_ua(2_IN,zetat) + &
                     & psi_ua(2_IN,zoq/obu) + &
                     & 0.8_r8*((-zetat)**(-onethird) - (-zeta)**(-onethird)) )
             else if (zeta.lt.0.0_r8) then
                ! Unstable regime
                ! EQN (12) with extra z0 term.
                re = loc_karman / &
                     & (log(zbot(n)/zoq) - psi_ua(2_IN,zeta) + psi_ua(2_IN,zoq/obu))
             else if (zeta.le.1.0_r8) then
                ! Stable regime
                ! EQN (13) with extra z0 term.
                re = loc_karman / &
                     & (log(zbot(n)/zoq) + 5.0_r8*zeta - 5.0_r8*zoq/obu)
             else
                ! Very stable regime
                ! EQN (14) with extra z0 term.
                re = loc_karman / &
                     & (log(obu/zoq) + 5.0_r8 - 5.0_r8*zoq/obu + &
                     & (5.0_r8*log(zeta) + zeta - 1.0_r8) )
             endif
             qstar = re * delq

             ! Update Obukhov length.
             thvstar = tstar*(1.0_r8 + 0.61_r8*qbot(n)) + 0.61_r8*thbot(n)*qstar
             ! EQN (4)
             obu = ustar*ustar * thv / (loc_karman*loc_g*thvstar)
             obu = sign( max(zbot(n)/10.0_r8, abs(obu)) ,obu)

             ! Update wind speed if in unstable regime.
             if (delthv.lt.0.0_r8) then
                ! EQN (20)
                wstar = beta * (-loc_g*ustar*thvstar*zpbl/thv)**onethird
                ! EQN (19)
                vmag = sqrt(vmag_rel**2 + wstar*wstar)
             else
                ! EQN (18)
                vmag = max(umin,vmag_rel)
             endif

          enddo ! End of iterations for ustar, tstar, qstar etc.


          !-----Calculate fluxes and wind stress.---------------------

          !--- momentum flux ---
          ! This should ensure zero wind stress when (relative) wind speed is zero,
          ! components are consistent with total, and we don't ever divide by zero.
          ! EQN (21)
          tau = rbot(n) * ustar * ustar
          taux(n) = tau * (ubot(n)-us(n)) / max(umin, vmag_rel)
          tauy(n) = tau * (vbot(n)-vs(n)) / max(umin, vmag_rel)

          !--- heat flux ---
          ! EQNs (22) and (23)
          sen (n) =  cp * rbot(n) * tstar * ustar
          lat (n) = xlv * rbot(n) * qstar * ustar
          lwup(n) = -loc_stebol * ts(n)**4

          !--- water flux ---
          evap(n) = lat(n)/xlv

          !---water isotope flux ---
          call wiso_flxoce(2,rbot(n),zbot(n),s16O(n),ts(n),r16O(n),ustar,re,ssq,evap_16O(n), &
               qbot(n),evap(n))
          call wiso_flxoce(3,rbot(n),zbot(n),sHDO(n),ts(n),rHDO(n),ustar,re,ssq, evap_HDO(n),&
               qbot(n),evap(n))
          call wiso_flxoce(4,rbot(n),zbot(n),s18O(n),ts(n),r18O(n),ustar,re,ssq, evap_18O(n), &
               qbot(n),evap(n))

          !------------------------------------------------------------
          ! compute diagnositcs: 2m ref T & Q, 10m wind speed squared
          !------------------------------------------------------------

          zeta = zbot(n) / obu
          if (zeta.lt.zetat) then
             if (zeta.lt.zetam) then
                ! Very unstable regime for U.
                ! EQN (7)
                u10 = vmag_abs + (ustar/loc_karman) * &
                     & 1.14_r8 * ((-zref/obu)**onethird - (-zeta)**onethird)
             else
                ! Unstable regime for U.
                ! EQN (8)
                u10 = vmag_abs + (ustar/loc_karman) * &
                     & (log(zref/zbot(n)) - (psi_ua(1_IN,zref/obu) - psi_ua(1_IN,zeta)) )
             endif
             ! Very unstable regime for T and q.
             ! EQN (11)
             tref(n) = thbot(n) + (tstar/loc_karman) * &
                  & 0.8_r8 * ((-zeta)**(-onethird) - (-ztref/obu)**(-onethird))
             qref(n) = qbot(n) + (qstar/loc_karman) * &
                  & 0.8_r8 * ((-zeta)**(-onethird) - (-ztref/obu)**(-onethird))

          else if (zeta.lt.0.0_r8) then
             ! Unstable regime.
             ! EQN (8)
             u10 = vmag_abs + (ustar/loc_karman) * &
                  & (log(zref/zbot(n)) - (psi_ua(1_IN,zref/obu) - psi_ua(1_IN,zeta)) )
             ! EQN (12)
             tref(n) = thbot(n) + (tstar/loc_karman) * &
                  & (log(ztref/zbot(n)) - (psi_ua(2_IN,ztref/obu) - psi_ua(2_IN,zeta)) )
             qref(n) = qbot(n) + (qstar/loc_karman) * &
                  & (log(ztref/zbot(n)) - (psi_ua(2_IN,ztref/obu) - psi_ua(2_IN,zeta)) )
          else if (zeta.le.1.0_r8) then
             ! Stable regime.
             ! EQN (9)
             u10 = vmag_abs + (ustar/loc_karman) * &
                  & (log(zref/zbot(n)) + 5.0_r8*zref/obu - 5.0_r8*zeta)
             ! EQN (13)
             tref(n) = thbot(n) + (tstar/loc_karman) * &
                  & (log(ztref/zbot(n)) + 5.0_r8*ztref/obu - 5.0_r8*zeta)
             qref(n) = qbot(n) + (qstar/loc_karman) * &
                  & (log(ztref/zbot(n)) + 5.0_r8*ztref/obu - 5.0_r8*zeta)
          else
             ! Very stable regime.
             ! EQN (10)
             u10 = vmag_abs + (ustar/loc_karman) * &
                  & (5.0_r8*log(zref/zbot(n)) + zref/obu - zeta)
             ! EQN (14)
             tref(n) = thbot(n) + (tstar/loc_karman) * &
                  & (5.0_r8*log(ztref/zbot(n)) + ztref/obu - zeta)
             qref(n) = qbot(n) + (qstar/loc_karman) * &
                  & (5.0_r8*log(ztref/zbot(n)) + ztref/obu - zeta)

          endif

          tref(n) = tref(n) - gamma*ztref   ! pot. temp to temp correction
          duu10n(n) = u10*u10 ! 10m wind speed squared

          !------------------------------------------------------------
          ! optional diagnostics, needed for water tracer fluxes (dcn)
          !------------------------------------------------------------
          if (present(ustar_sv)) ustar_sv(n) = ustar
          if (present(ssq_sv  )) ssq_sv(n)   = ssq
          if (present(re_sv   )) re_sv(n)    = re


       else

          !------------------------------------------------------------
          ! no valid data here -- out of ocean domain
          !------------------------------------------------------------
          sen   (n) = spval  ! sensible         heat flux  (W/m^2)
          lat   (n) = spval  ! latent           heat flux  (W/m^2)
          lwup  (n) = spval  ! long-wave upward heat flux  (W/m^2)
          evap  (n) = spval  ! evaporative water flux ((kg/s)/m^2)
          evap_16O (n) = spval !water tracer flux (kg/s)/m^2)
          evap_HDO (n) = spval !HDO tracer flux  (kg/s)/m^2)
          evap_18O (n) = spval !H218O tracer flux (kg/s)/m^2)
          taux  (n) = spval  ! x surface stress (N)
          tauy  (n) = spval  ! y surface stress (N)
          tref  (n) = spval  !  2m reference height temperature (K)
          qref  (n) = spval  !  2m reference height humidity (kg/kg)
          duu10n(n) = spval  ! 10m wind speed squared (m/s)^2
          ! Optional diagnostics too:
          if (present(ustar_sv)) ustar_sv(n) = spval
          if (present(re_sv   )) re_sv   (n) = spval
          if (present(ssq_sv  )) ssq_sv  (n) = spval

       endif

    ENDDO ! loop over grid points

  END subroutine flux_atmOcn_UA

  !===============================================================================
  ! Functions/subroutines used by UA surface flux scheme.
  !===============================================================================

  ! Stability function for rb < 0

  real(r8) function psi_ua(k,zeta)

    implicit none

    !-----Input variables.----------
    integer(in), intent(in) :: k       ! Indicates whether this is for momentum (k=1)
    ! or for heat/moisture (k=2)
    real(r8), intent(in) :: zeta       ! Dimensionless height (=z/L)

    !-----Local variables.----------
    real(r8) :: chik                   ! Function of zeta.

    ! EQN (16)
    chik = (1.0_r8 - 16.0_r8*zeta)**0.25_r8

    if(k.eq.1) then
       ! EQN (15) for momentum
       psi_ua = 2.0_r8 * log((1.0_r8 + chik)*0.5_r8) + &
            &      log((1.0_r8 + chik*chik)*0.5_r8) - &
            & 2.0_r8 * atan(chik) + 2.0_r8 * atan(1.0_r8)
    else
       ! EQN (15) for heat/moisture
       psi_ua = 2.0_r8 * log((1.0_r8 + chik*chik)*0.5_r8)
    endif

  end function psi_ua

  !===============================================================================
  ! Uses Tetens' formula for saturation vapor pressure from
  ! Buck(1981) JAM 20, 1527-1532

  real(r8) function qsat_ua(t,p,loc_epsilon)

    implicit none

    !-----Input variables.----------
    real(r8), intent(in) :: t           ! temperature (K)
    real(r8), intent(in) :: p           ! pressure (Pa)
    real(r8), intent(in) :: loc_epsilon ! Ratio of gas constants (-)

    !-----Local variables.----------
    real(r8) :: esat                    ! saturated vapor pressure (hPa)

    ! Calculate saturated vapor pressure in hPa.
    esat = (1.0007_r8 + 0.00000346_r8 * (p/100.0_r8)) * 6.1121_r8 * &
         & exp(17.502_r8 * (t - loc_tkfrz) / (240.97_r8 + (t - loc_tkfrz)))

    ! Convert to specific humidity (kg kg-1).
    qsat_ua = loc_epsilon * esat / ((p/100.0_r8) - (1.0_r8 - loc_epsilon)*esat)

  end function qsat_ua

  !===============================================================================
  ! Calculate roughness lengths: zo, zot, zoq.

  subroutine rough_ua(zo,zot,zoq,ustar,visa)

    implicit none

    !-----Input variables.----------
    real(r8), intent(in) :: ustar      ! friction velocity (m s-1)
    real(r8), intent(in) :: visa       ! kinematic viscosity of dry air (m2 s-1)

    !-----Output variables.---------
    real(r8), intent(out) :: zo        ! roughness length for momentum (m)
    real(r8), intent(out) :: zot       ! roughness length for heat (m)
    real(r8), intent(out) :: zoq       ! roughness length for water vapor (m)

    !-----Local variables.----------
    real(r8) :: re_rough               ! Rougness Reynold's number (-)
    real(r8) :: xq                     ! Logarithm of roughness length ratios (moisture)
    real(r8) :: xt                     ! Logarithm of roughness length ratios (heat)

    zo = 0.013_r8*ustar*ustar/loc_g + 0.11_r8*visa/ustar      ! EQN (24)
    re_rough = ustar*zo/visa                                  ! By definition.
    xq = 2.67_r8*re_rough**0.25_r8 - 2.57_r8                  ! EQN (25)
    xt = xq                                                   ! EQN (26)
    zoq = zo/exp(xq)                                          ! By definition of xq
    zot = zo/exp(xt)                                          ! By definition of xt

  end subroutine rough_ua

  real(r8) elemental function cuberoot(a)
    real(r8), intent(in) :: a
    real(r8), parameter :: one_third = 1._r8/3._r8
    cuberoot = sign(abs(a)**one_third, a)
  end function cuberoot

  !===============================================================================
  ! !IROUTINE: flux_atmOcn_diurnal -- internal atm/ocn flux calculation
  !
  ! !DESCRIPTION:
  !
  !     Internal atm/ocn flux calculation
  !
  ! !REVISION HISTORY:
  !     2002-Jun-10 - B. Kauffman - code migrated from cpl5 to cpl6
  !     2003-Apr-02 - B. Kauffman - taux & tauy now utilize ocn velocity
  !     2003-Apr-02 - B. Kauffman - tref,qref,duu10n mods as per Bill Large
  !     2006-Nov-07 - B. Kauffman - code migrated from cpl6 to share
  !===============================================================================
  SUBROUTINE flux_atmOcn_diurnal &
       (logunit, nMax  ,zbot  ,ubot  ,vbot  ,thbot ,   &
       qbot  ,s16O  ,sHDO  ,s18O  ,rbot  ,             &
       tbot  ,us    ,vs    ,                           &
       ts    ,mask  , seq_flux_atmocn_minwind,         &
       sen   ,lat   ,lwup  ,                           &
       r16O  ,rhdo  ,r18O  ,evap  ,evap_16O,           &
       evap_HDO     ,evap_18O,                         &
       taux  ,tauy  ,tref  ,qref  ,                    &
       uGust, lwdn , swdn , swup, prec   ,             &
       swpen, ocnsal, ocn_prognostic, flux_diurnal,    &
       ocn_surface_flux_scheme,                                    &
       latt, long , warm , salt , speed, regime,       &
       warmMax, windMax, qSolAvg, windAvg,             &
       warmMaxInc, windMaxInc, qSolInc, windInc, nInc, &
       tBulk, tSkin, tSkin_day, tSkin_night,           &
       cSkin, cSkin_night, secs ,dt,                   &
       duu10n,  ustar_sv   ,re_sv ,ssq_sv,             &
       missval, cold_start    )
    ! !USES:

    use water_isotopes, only: wiso_flxoce !subroutine used to calculate water isotope fluxes.

    implicit none

    ! !INPUT/OUTPUT PARAMETERS:

    !--- input arguments --------------------------------
    integer    ,intent(in) :: logunit
    integer(in),intent(in) :: nMax        ! data vector length
    integer(in),intent(in) :: mask (nMax) ! ocn domain mask       0 <=> out of domain
    real(r8)   ,intent(in) :: zbot (nMax) ! atm level height      (m)
    real(r8)   ,intent(in) :: ubot (nMax) ! atm u wind            (m/s)
    real(r8)   ,intent(in) :: vbot (nMax) ! atm v wind            (m/s)
    real(r8)   ,intent(in) :: thbot(nMax) ! atm potential T       (K)
    real(r8)   ,intent(in) :: qbot (nMax) ! atm specific humidity (kg/kg)
    real(r8)   ,intent(in) :: s16O (nMax) ! atm H216O tracer conc. (kg/kg)
    real(r8)   ,intent(in) :: sHDO (nMax) ! atm HDO tracer conc.  (kg/kg)
    real(r8)   ,intent(in) :: s18O (nMax) ! atm H218O tracer conc. (kg/kg)
    real(r8)   ,intent(in) :: r16O (nMax) ! ocn H216O tracer ratio/Rstd
    real(r8)   ,intent(in) :: rHDO (nMax) ! ocn HDO tracer ratio/Rstd
    real(r8)   ,intent(in) :: r18O (nMax) ! ocn H218O tracer ratio/Rstd
    real(r8)   ,intent(in) :: rbot (nMax) ! atm air density       (kg/m^3)
    real(r8)   ,intent(in) :: tbot (nMax) ! atm T                 (K)
    real(r8)   ,intent(in) :: us   (nMax) ! ocn u-velocity        (m/s)
    real(r8)   ,intent(in) :: vs   (nMax) ! ocn v-velocity        (m/s)
    real(r8)   ,intent(in) :: ts   (nMax) ! ocn temperature       (K)

    !--- new    arguments -------------------------------
    real(r8),intent(inout) :: swpen (nMax)       ! NEW
    real(r8),intent(inout) :: ocnsal(nMax)       ! NEW (kg/kg)
    logical ,intent(in)    :: ocn_prognostic     ! NEW
    logical ,intent(in)    :: flux_diurnal       ! NEW logical for diurnal on/off
    integer(in) ,intent(in)    :: ocn_surface_flux_scheme

    real(r8),intent(in)    :: uGust (nMax)      ! NEW not used
    real(r8),intent(in)    :: lwdn  (nMax)       ! NEW
    real(r8),intent(in)    :: swdn  (nMax)       ! NEW
    real(r8),intent(in)    :: swup  (nMax)       ! NEW
    real(r8),intent(in)    :: prec  (nMax)       ! NEW
    real(r8),intent(in)    :: latt  (nMax)       ! NEW
    real(r8),intent(in)    :: long  (nMax)       ! NEW
    real(r8),intent(inout) :: warm  (nMax)       ! NEW
    real(r8),intent(inout) :: salt  (nMax)       ! NEW
    real(r8),intent(inout) :: speed (nMax)       ! NEW
    real(r8),intent(inout) :: regime(nMax)       ! NEW
    real(r8),intent(out)   :: warmMax(nMax)      ! NEW
    real(r8),intent(out)   :: windMax(nMax)      ! NEW
    real(r8),intent(inout) :: qSolAvg(nMax)      ! NEW
    real(r8),intent(inout) :: windAvg(nMax)      ! NEW
    real(r8),intent(inout) :: warmMaxInc(nMax)   ! NEW
    real(r8),intent(inout) :: windMaxInc(nMax)   ! NEW
    real(r8),intent(inout) :: qSolInc(nMax)      ! NEW
    real(r8),intent(inout) :: windInc(nMax)      ! NEW
    real(r8),intent(inout) :: nInc(nMax)         ! NEW

    real(r8),intent(out)   :: tBulk (nMax)       ! NEW
    real(r8),intent(out)   :: tSkin (nMax)       ! NEW
    real(r8),intent(out)   :: tSkin_day (nMax)   ! NEW
    real(r8),intent(out)   :: tSkin_night (nMax) ! NEW
    real(r8),intent(out)   :: cSkin (nMax)       ! NEW
    real(r8),intent(out)   :: cSkin_night (nMax) ! NEW
    integer(in),intent(in) :: secs               ! NEW  elsapsed seconds in day (GMT)
    integer(in),intent(in) :: dt                 ! NEW
    logical ,intent(in)    :: cold_start         ! cold start flag
    real(r8),intent(in)    :: seq_flux_atmocn_minwind   ! minimum wind speed for atmocn      (m/s)

    real(r8),intent(in) ,optional :: missval     ! masked value

    !--- output arguments -------------------------------
    real(r8),intent(out)  ::  sen  (nMax) ! heat flux: sensible    (W/m^2)
    real(r8),intent(out)  ::  lat  (nMax) ! heat flux: latent      (W/m^2)
    real(r8),intent(out)  ::  lwup (nMax) ! heat flux: lw upward   (W/m^2)
    real(r8),intent(out)  ::  evap (nMax) ! water flux: evap  ((kg/s)/m^2)
    real(r8),intent(out)  ::  evap_16O (nMax) ! water flux: evap ((kg/s/m^2)
    real(r8),intent(out)  ::  evap_HDO (nMax) ! water flux: evap ((kg/s)/m^2)
    real(r8),intent(out)  ::  evap_18O (nMax) ! water flux: evap ((kg/s/m^2)
    real(r8),intent(out)  ::  taux (nMax) ! surface stress, zonal      (N)
    real(r8),intent(out)  ::  tauy (nMax) ! surface stress, maridional (N)
    real(r8),intent(out)  ::  tref (nMax) ! diag:  2m ref height T     (K)
    real(r8),intent(out)  ::  qref (nMax) ! diag:  2m ref humidity (kg/kg)
    real(r8),intent(out)  :: duu10n(nMax) ! diag: 10m wind speed squared (m/s)^2

    real(r8),intent(out),optional :: ustar_sv(nMax) ! diag: ustar
    real(r8),intent(out),optional :: re_sv   (nMax) ! diag: sqrt of exchange coefficient (water)
    real(r8),intent(out),optional :: ssq_sv  (nMax) ! diag: sea surface humidity  (kg/kg)


    !--- local constants --------------------------------
    real(r8),parameter :: zref  = 10.0_r8 ! reference height           (m)
    real(r8),parameter :: ztref =  2.0_r8 ! reference height for air T (m)

    real(r8),parameter :: lambdaC  = 6.0_r8
    real(r8),parameter :: lambdaL  = 0.0_r8
    real(r8),parameter :: doLMax   = 1.0_r8
    real(r8),parameter :: pwr      = 0.2_r8
    real(r8),parameter :: Rizero   = 1.0_r8
    real(r8),parameter :: NUzero   = 40.0e-4_r8
    real(r8),parameter :: Prandtl  = 1.0_r8
    real(r8),parameter :: kappa0   = 0.2e-4_r8

    real(r8),parameter :: F0       = 0.5_r8
    real(r8),parameter :: F1       = 0.15_r8
    real(r8),parameter :: R1       = 10.0_r8

    real(r8),parameter :: Ricr     = 0.30_r8
    real(r8),parameter :: tiny     = 1.0e-12_r8
    real(r8),parameter :: tiny2    = 1.0e-6_r8
    real(r8),parameter :: pi       = SHR_CONST_PI

    !!++ COARE only
    real(r8),parameter :: zpbl =700.0_r8 ! PBL depth [m] for gustiness parametriz.

    !--- local variables --------------------------------
    integer(in) :: n       ! vector loop index
    integer(in) :: iter       ! iteration loop index
    integer(in) :: lsecs   ! local seconds elapsed
    integer(in) :: lonsecs ! incrememnt due to lon offset
    real(r8)    :: vmag    ! surface wind magnitude   (m/s)
    real(r8)    :: ssq     ! sea surface humidity     (kg/kg)
    real(r8)    :: delt    ! potential T difference   (K)
    real(r8)    :: delq    ! humidity difference      (kg/kg)
    real(r8)    :: stable  ! stability factor
    real(r8)    :: rdn     ! sqrt of neutral exchange coeff (momentum)
    real(r8)    :: rhn     ! sqrt of neutral exchange coeff (heat)
    real(r8)    :: ren     ! sqrt of neutral exchange coeff (water)
    real(r8)    :: rd      ! sqrt of exchange coefficient (momentum)
    real(r8)    :: rh      ! sqrt of exchange coefficient (heat)
    real(r8)    :: re      ! sqrt of exchange coefficient (water)
    real(r8)    :: ustar   ! ustar
    real(r8)    :: ustar_prev   ! ustar
    real(r8)    :: qstar   ! qstar
    real(r8)    :: tstar   ! tstar
    real(r8)    :: hol     ! H (at zbot) over L
    real(r8)    :: xsq     ! ?
    real(r8)    :: xqq     ! ?
    real(r8)    :: psimh   ! stability function at zbot (momentum)
    real(r8)    :: psixh   ! stability function at zbot (heat and water)
    real(r8)    :: psix2   ! stability function at ztref reference height
    real(r8)    :: alz     ! ln(zbot/zref)
    real(r8)    :: al2     ! ln(zref/ztref)
    real(r8)    :: u10n    ! 10m neutral wind
    real(r8)    :: tau     ! stress at zbot
    real(r8)    :: cp      ! specific heat of moist air
    real(r8)    :: fac     ! vertical interpolation factor
    real(r8)    :: DTiter  !
    real(r8)    :: DSiter  !
    real(r8)    :: DViter  !

    real(r8)    :: Dcool   !
    real(r8)    :: Qdel    ! net cool skin heating
    real(r8)    :: Hd      ! net heating above -z=d
    real(r8)    :: Hb      ! net kinematic heating above -z = delta
    real(r8)    :: lambdaV !
    real(r8)    :: Fd      ! net fresh water forcing above -z=d
    real(r8)    :: ustarw  ! surface wind forcing of layer above -z=d

    real(r8)    :: Qsol   ! solar heat flux (W/m2)
    real(r8)    :: Qnsol  ! non-solar heat flux (W/m2)

    real(r8)    :: SSS  ! sea surface salinity
    real(r8)    :: alphaT  !
    real(r8)    :: betaS  !

    real(r8)    :: doL     ! ocean forcing stablity parameter
    real(r8)    :: Rid     ! Richardson number at depth d
    real(r8)    :: Ribulk  ! Bulk  Richardson number at depth d
    real(r8)    :: FofRi   ! Richardon number dependent diffusivity
    real(r8)    :: Smult   ! multiplicative term based on regime
    real(r8)    :: Sfact   ! multiplicative term based on regime
    real(r8)    :: Kdiff   ! diffusive term based on regime
    real(r8)    :: Kvisc   ! viscosity term based on regime
    real(r8)    :: rhocn   !
    real(r8)    :: rcpocn  !
    real(r8)    :: Nreset  ! value for multiplicative reset factor
    logical     :: lmidnight
    logical     :: ltwopm
    logical     :: ltwoam
    logical     :: lfullday
    integer     :: nsum
    real(r8)    :: pexp   ! eqn 19
    real(r8)    :: AMP    ! eqn 18
    real(r8)    :: dif3
    real(r8)    :: phid
    real(r8)    :: spval

    !!++ COARE only
    real(r8)    :: zo,zot,zoq      ! roughness lengths
    real(r8)    :: hsb,hlb         ! sens & lat heat flxs at zbot
    real(r8)    :: trf,qrf,urf,vrf ! reference-height quantities

    !--- local functions --------------------------------
    real(r8)    :: qsat   ! function: the saturation humididty of air (kg/m^3)
    real(r8)    :: cdn    ! function: neutral drag coeff at 10m
    real(r8)    :: psimhu ! function: unstable part of psimh
    real(r8)    :: psixhu ! function: unstable part of psimx
    real(r8)    :: Umps   ! dummy arg ~ wind velocity (m/s)
    real(r8)    :: Tk     ! dummy arg ~ temperature (K)
    real(r8)    :: xd     ! dummy arg ~ ?
    real(r8)    :: molvisc ! molecular viscosity
    real(r8)    :: molPr   ! molecular Prandtl number

    !--- for cold air outbreak calc --------------------------------
    real(r8)    :: tdiff(nMax)               ! tbot - ts
    real(r8)    :: vscl

    qsat(Tk)   = 640380.0_r8 / exp(5107.4_r8/Tk)
    cdn(Umps)  =   0.0027_r8 / Umps + 0.000142_r8 + 0.0000764_r8 * Umps
    psimhu(xd) = log((1.0_r8+xd*(2.0_r8+xd))*(1.0_r8+xd*xd)/8.0_r8) - 2.0_r8*atan(xd) + 1.571_r8
    psixhu(xd) = 2.0_r8 * log((1.0_r8 + xd*xd)/2.0_r8)
    molvisc(Tk)  = 1.623e-6_r8 * exp((-1.0_r8*(Tk-273.15_r8))/45.2_r8)
    molPr(Tk)    = 11.64_r8 * exp((-1.0_r8*(Tk-273.15_r8))/40.7_r8)

    !--- formats ----------------------------------------
    character(*),parameter :: subName = '(flux_atmOcn_diurnal) '
    character(*),parameter ::   F00 = "('(flux_atmOcn_diurnal) ',4a)"

    !-------------------------------------------------------------------------------
    ! PURPOSE:
    !   computes atm/ocn surface fluxes
    !
    ! NOTES:
    !   o all fluxes are positive downward
    !   o net heat flux = net sw + lw up + lw down + sen + lat
    !   o here, tstar = <WT>/U*, and qstar = <WQ>/U*.
    !   o wind speeds should all be above a minimum speed (eg. 1.0 m/s)
    !
    ! ASSUMPTIONS:
    !   o Neutral 10m drag coeff: cdn = .0027/U10 + .000142 + .0000764 U10
    !   o Neutral 10m stanton number: ctn = .0327 sqrt(cdn), unstable
    !                                 ctn = .0180 sqrt(cdn), stable
    !   o Neutral 10m dalton number:  cen = .0346 sqrt(cdn)
    !   o The saturation humidity of air at T(K): qsat(T)  (kg/m^3)
    !-------------------------------------------------------------------------------

    if (debug > 0) write(logunit,F00) "enter"

    ! this is especially for flux_diurnal calculations
    if (.not. flux_diurnal) then
       write(logunit,F00) "ERROR: flux_diurnal must be true"
       call shr_sys_abort(subName//"flux diurnal must be true")
    endif
    spval = shr_const_spval
    rh = spval
    dviter = spval
    dtiter = spval
    dsiter = spval
    al2 = log(zref/ztref)
    !--- for cold air outbreak calc --------------------------------
    tdiff= tbot - ts

    ! equations 18 and 19
    AMP = 1.0_r8/F0-1.0_r8
    pexp = log( (1.0_r8/F1-F0) / (1.0_r8-F0) ) / log(R1)

    if (.not. ocn_prognostic) then
       ! Set swpen and ocean salinity from following analytic expressions
       swpen(:) = 0.67_r8*(exp((-1._r8*shr_const_zsrflyr)/1.0_r8)) + &
            0.33_r8*exp((-1._r8*shr_const_zsrflyr)/17.0_r8)
       ocnsal(:) = shr_const_ocn_ref_sal/1000.0_r8
    else
       ! use swpen and ocnsal from input argument
    endif

    if (cold_start) then
       write(logunit,F00) "Initialize diurnal cycle fields"
       warm       (:) = 0.0_r8
       salt       (:) = 0.0_r8
       speed      (:) = 0.0_r8
       regime     (:) = 0.0_r8
       qSolAvg    (:) = 0.0_r8
       windAvg    (:) = 0.0_r8
       warmMax    (:) = 0.0_r8
       windMax    (:) = 0.0_r8
       warmMaxInc (:) = 0.0_r8
       windMaxInc (:) = 0.0_r8
       qSolInc    (:) = 0.0_r8
       windInc    (:) = 0.0_r8
       nInc       (:) = 0.0_r8
       tSkin_day  (:) = ts(:)
       tSkin_night(:) = ts(:)
       cSkin_night(:) = 0.0_r8
    endif
    u10n = 0.0_r8
    stable = 0.0_r8
    DO n=1,nMax

       if (mask(n) /= 0) then

          !--- compute some initial and useful flux quantities ---

          vmag     = max(seq_flux_atmocn_minwind, sqrt( (ubot(n)-us(n))**2 + (vbot(n)-vs(n))**2) )
          if (use_coldair_outbreak_mod) then
             ! Cold Air Outbreak Modification:
             ! Increase windspeed for negative tbot-ts
             ! based on Mahrt & Sun 1995,MWR

             if (tdiff(n).lt.td0) then
                vscl=min((1._r8+alpha*(abs(tdiff(n)-td0)**0.5_r8/abs(vmag))),maxscl)
                vmag=vmag*vscl
             endif
          endif
          alz      = log(zbot(n)/zref)
          hol      = 0.0
          psimh    = 0.0
          psixh    = 0.0
          rdn      = sqrt(cdn(vmag))

          tBulk(n) = ts(n)+warm(n)    ! first guess for tBulk from read in ts,warm
          tSkin(n) = tBulk(n)
          Qsol     = swdn(n) + swup(n)
          SSS      = 1000.0_r8*ocnsal(n)+salt(n)
          lambdaV = lambdaC

          alphaT   = 0.000297_r8*(1.0_r8+0.0256_r8*(ts(n)-298.15_r8)+0.003_r8*(SSS - 35.0_r8))
          betaS    = 0.000756_r8*(1.0_r8-0.0016_r8*(ts(n)-298.15_r8))
          rhocn    = 1023.342_r8*(1.0_r8-0.000297_r8*(ts(n)-298.15_r8)+0.000756_r8 * (SSS - 35.0_r8))
          rcpocn   = rhocn * 3990.0_r8*(1.0_r8-0.0012_r8*(SSS - 35.0_r8))

          Rid =  shr_const_g * (alphaT*warm(n) - betaS*salt(n)) *pwr*shr_const_zsrflyr  / &
               ( pwr*MAX(tiny,speed(n)) )**2

          Ribulk = 0.0

          !----------------------------------------------------------
          ! convert elapsed time from GMT to local &
          ! check elapsed time. reset warm if near lsecs = reset_sec
          !----------------------------------------------------------
          Nreset = 1.0_r8

          lonsecs   = ceiling(long(n)/360.0_r8*86400.0)
          lsecs     = mod(secs + lonsecs,86400)

          lmidnight = (lsecs >= 0     .and. lsecs < dt)        ! 0 = midnight
          ltwopm    = (lsecs >= 48600 .and. lsecs < 48600+dt)  ! 48600 = 1:30pm
          ltwoam    = (lsecs >= 5400  .and. lsecs < 5400 +dt)  ! 5400 = 1:30am
          lfullday  = (lsecs > 86400-dt .and. lsecs <= 86400)
          nsum = nint(nInc(n))

          if ( lmidnight ) then
             Regime(n)  = 1.0_r8               !  RESET DIURNAL
             warm(n)    = 0.0_r8
             salt(n)    = 0.0_r8
             speed(n)   = 0.0_r8
          endif

          ssq    = 0.98_r8 * qsat(tBulk(n)) / rbot(n)   ! sea surf hum (kg/kg)
          delt   = thbot(n) - tBulk(n)                  ! pot temp diff (K)
          delq   = qbot(n) - ssq                     ! spec hum dif (kg/kg)
          cp     = shr_const_cpdair*(1.0_r8 + shr_const_cpvir*ssq)

          !!.................................................................
          !! ocn_surface_flux_scheme = 0 : Default E3SMv1
          !!                         = 1 : COARE algorithm
          !!.................................................................
          if (ocn_surface_flux_scheme .eq. 0) then! use Large algorithm
             stable = 0.5_r8 + sign(0.5_r8 , delt)


             !--- shift wind speed using old coefficient  and stability function

             rd   = rdn / (1.0_r8 + rdn/shr_const_karman*(alz-psimh))
             u10n = vmag * rd / rdn

             !--- initial neutral  transfer coeffs at 10m
             rdn    = sqrt(cdn(u10n))
             rhn    = (1.0_r8-stable) * 0.0327_r8 + stable * 0.018_r8
             ren    = 0.0346_r8

             !--- initial ustar, tstar, qstar ---
             ustar = rdn * vmag
             tstar = rhn * delt
             qstar = ren * delq

          else if (ocn_surface_flux_scheme .eq. 1) then! use COARE algorithm

             call cor30a(ubot(n),vbot(n),tbot(n),qbot(n),rbot(n) &  ! in atm params
                  & ,us(n),vs(n),tBulk(n),ssq                &  ! in surf params (NB ts -> tBulk)
                  & ,zpbl,zbot(n),zbot(n),zref,ztref,ztref   &  ! in heights
                  & ,tau,hsb,hlb                             &  ! out: fluxes
                  & ,zo,zot,zoq,hol,ustar,tstar,qstar        &  ! out: ss scales
                  & ,rd,rh,re                                &  ! out: exch. coeffs
                  & ,trf,qrf,urf,vrf)                             ! out: reference-height params
             ! for the sake of maintaining same defs
             hol=zbot(n)/hol
             rd=sqrt(rd)
             rh=sqrt(rh)
             re=sqrt(re)

          ELSE  ! N.B.: *no* valid ocn_surface_flux_scheme=2 option if diurnal=.true.

             call shr_sys_abort(subName//" flux_atmOcn_diurnal requires ocn_surface_flux_scheme = 0 or 1")
          ENDIF

          ustar_prev = ustar * 2.0_r8
          iter = 0
          ! --- iterate ---
          ! Originally this code did three iterations while the non-diurnal version did two
          ! So in the new loop this is <= flux_con_max_iter instead of < so that the same defaults
          ! will give the same answers in both cases.
          do while( abs((ustar - ustar_prev)/ustar) > flux_con_tol .and. iter <= flux_con_max_iter)
             iter = iter + 1
             ustar_prev = ustar
             !------------------------------------------------------------
             ! iterate to converge on FLUXES  Z/L, ustar, tstar and qstar
             ! and on Rid  in the DIURNAL CYCLE
             !------------------------------------------------------------
             Smult = 0.0_r8
             Sfact = 0.0_r8
             Kdiff = 0.0_r8
             Kvisc = 0.0_r8
             dif3 = 0.0_r8

             ustarw  = ustar*sqrt(max(tiny,rbot(n)/rhocn))
             Qnsol   = lwdn(n) - shr_const_stebol*(tSkin(n))**4 + &
                  rbot(n)*ustar*(cp*tstar + shr_const_latvap*qstar)
             Hd      = (Qnsol   + Qsol*(1.0_r8-swpen(n)) ) / rcpocn
             Fd      = (prec(n) + rbot(n)*ustar*qstar ) * SSS / rhocn

             !--- COOL SKIN EFFECT ---
             Dcool  = lambdaV*molvisc(tBulk(n)) / ustarw
             Qdel   = Qnsol + Qsol * &
                  (0.137_r8 + 11.0_r8*Dcool - 6.6e-5/Dcool *(1.0_r8 - exp((-1.0_r8*Dcool)/8.0e-4)))
             Hb = (Qdel/rcpocn)+(Fd*betaS/alphaT)
             Hb = min(Hb , 0.0_r8)

             !            lambdaV = lambdaC*(1.0_r8 + ( (0.0_r8-Hb)*16.0_r8*molvisc(tBulk(n))* &
             !                 shr_const_g*alphaT*molPr(tBulk(n))**2/ustarw**4)**0.75)**(-1._r8/3._r8)
             lambdaV = 6.5_r8
             cSkin(n) =  MIN(0.0_r8, lambdaV * molPr(tBulk(n)) * Qdel / ustarw / rcpocn )

             !--- REGIME ---
             doL = shr_const_zsrflyr*shr_const_karman*shr_const_g* &
                  (alphaT*Hd + betaS*Fd ) / ustarw**3
             Rid = MAX(0.0_r8,Rid)
             Smult = dt * (pwr+1.0_r8) / (shr_const_zsrflyr*pwr)
             Sfact = dt * (pwr+1.0_r8) / (shr_const_zsrflyr)**2
             FofRi = 1.0_r8/(1.0_r8 + AMP*(Rid/Rizero)**pexp)

             if ( (doL.gt.0.0_r8) .and. (Qsol.gt.0.0)  ) then
                phid  = MIN(1.0_r8 + 5.0_r8 * doL, 5.0_r8 + doL)
                FofRi = 1.0_r8/(1.0_r8 + AMP*(Rid/Rizero)**pexp)
                dif3 = (kappa0 + NUzero *FofRi)

                if ((doL.le.lambdaL).and.(NINT(regime(n)).le.2)) then
                   regime(n) = 2.0_r8
                   Kdiff =  shr_const_karman * ustarw * shr_const_zsrflyr / phid
                   Kvisc = Kdiff * (1.0_r8 - doL/lambdaL)**2 + &
                        dif3 * (doL/lambdaL)**2 * (3.0_r8 - 2.0_r8 * doL/lambdaL)
                   Kdiff = Kvisc
                else
                   regime(n) = 3.0_r8
                   Kdiff =          kappa0 + NUzero * FofRi
                   Kvisc = Prandtl* kappa0 + NUzero * FofRi
                endif
             else
                if (regime(n).eq.1.0_r8) then
                   Smult      = 0.0_r8
                else
                   if (Ribulk .gt. Ricr) then
                      regime(n) = 3.0_r8
                      Kdiff =          kappa0 + NUzero * FofRi
                      Kvisc = Prandtl* kappa0 + NUzero * FofRi
                   else
                      regime(n) = 4.0_r8
                      Kdiff = shr_const_karman*ustarw*shr_const_zsrflyr *cuberoot(1.0_r8-7.0_r8*doL)
                      Kvisc = Kdiff
                   endif
                endif

             endif

             !--- IMPLICIT INTEGRATION ---

             DTiter = (warm(n)  +(Smult*Hd))               /(1.+ Sfact*Kdiff)
             DSiter = (salt(n)  -(Smult*Fd))               /(1.+ Sfact*Kdiff)
             DViter = (speed(n) +(Smult*ustarw*ustarw))    /(1.+ Sfact*Kvisc)
             DTiter = MAX( 0.0_r8, DTiter)
             DViter = MAX( 0.0_r8, DViter)

             Rid =(shr_const_g*(alphaT*DTiter-betaS*DSiter)*pwr*shr_const_zsrflyr)  / &
                  (pwr*MAX(tiny,DViter))**2
             Ribulk = Rid * pwr
             Ribulk = 0.0_r8
             tBulk(n) = ts(n) + DTiter
             tSkin(n) = tBulk(n) + cskin(n)

             !--need to update ssq,delt,delq as function of tBulk ----

             ssq    = 0.98_r8 * qsat(tBulk(n)) / rbot(n)   ! sea surf hum (kg/kg)
             delt   = thbot(n) - tBulk(n)                  ! pot temp diff (K)
             delq   = qbot(n) - ssq                        ! spec hum dif (kg/kg)

             !--- UPDATE FLUX ITERATION ---

             !!.................................................................
             !! ocn_surface_flux_scheme = 0 : Default CESM1.2
             !!                         = 1 : COARE algorithm
             !!.................................................................
             if (ocn_surface_flux_scheme .eq. 0) then! use Large algorithm

                !--- compute stability & evaluate all stability functions ---
                hol  = shr_const_karman*shr_const_g*zbot(n)*  &
                     (tstar/thbot(n)+qstar/(1.0_r8/shr_const_zvir+qbot(n)))/ustar**2
                hol  = sign( min(abs(hol),10.0_r8), hol )
                stable = 0.5_r8 + sign(0.5_r8 , hol)
                xsq    = max(sqrt(abs(1.0_r8 - 16.0_r8*hol)) , 1.0_r8)
                xqq    = sqrt(xsq)
                psimh  = -5.0_r8*hol*stable + (1.0_r8-stable)*psimhu(xqq)
                psixh  = -5.0_r8*hol*stable + (1.0_r8-stable)*psixhu(xqq)

                !--- shift wind speed using old coefficient  and stability function  ---
                rd   = rdn / (1.0_r8 + rdn/shr_const_karman*(alz-psimh))
                u10n = vmag * rd / rdn

                !--- update neutral  transfer coeffs at 10m
                rdn    = sqrt(cdn(u10n))
                rhn    = (1.0_r8-stable) * 0.0327_r8 + stable * 0.018_r8
                ren    = 0.0346_r8

                !--- shift all coeffs to measurement height and stability ---
                rd = rdn / (1.0_r8 + rdn/shr_const_karman*(alz-psimh))
                rh = rhn / (1.0_r8 + rhn/shr_const_karman*(alz-psixh))
                re = ren / (1.0_r8 + ren/shr_const_karman*(alz-psixh))

                ustar = rd * vmag
                tstar = rh * delt
                qstar = re * delq

                !--- heat flux ---

                tau     = rbot(n) * ustar * ustar
                sen (n) =                cp * tau * tstar / ustar
                lat (n) = shr_const_latvap * tau * qstar / ustar

             else if (ocn_surface_flux_scheme .eq. 1) then! use COARE algorithm

                call cor30a(ubot(n),vbot(n),tbot(n),qbot(n),rbot(n) &  ! in atm params
                     & ,us(n),vs(n),tBulk(n),ssq                &  ! in surf params (NB ts -> tBulk)
                     & ,zpbl,zbot(n),zbot(n),zref,ztref,ztref   &  ! in heights
                     & ,tau,hsb,hlb                             &  ! out: fluxes
                     & ,zo,zot,zoq,hol,ustar,tstar,qstar        &  ! out: ss scales
                     & ,rd,rh,re                                &  ! out: exch. coeffs
                     & ,trf,qrf,urf,vrf)                               ! out: reference-height params
                ! for the sake of maintaining same defs
                hol=zbot(n)/hol
                rd=sqrt(rd)
                rh=sqrt(rh)
                re=sqrt(re)

                !--- heat flux ---

                sen (n) =  hsb
                lat (n) =  hlb

             else ! N.B.: NO ocn_surface_flux_scheme=2 option
                call shr_sys_abort(subName//", flux_diurnal requires ocn_surface_flux_scheme = 0 or 1")
             endif

          ENDDO   ! end iteration loop
          if (iter < 1) then
             call shr_sys_abort('No iterations performed ')
          end if
          !--- COMPUTE FLUXES TO ATMOSPHERE AND OCEAN ---

          ! Now calculated further up in subroutine.
          !tau = rbot(n) * ustar * ustar
          !sen (n) =                cp * tau * tstar / ustar
          !lat (n) =  shr_const_latvap * tau * qstar / ustar

          !--- momentum flux ---
          taux(n) = tau * (ubot(n)-us(n)) / vmag
          tauy(n) = tau * (vbot(n)-vs(n)) / vmag

          !--- LW radiation ---
          lwup(n) = -shr_const_stebol * Tskin(n)**4

          !--- water flux ---
          evap(n) = lat(n)/shr_const_latvap

          !---water isotope flux ---
          !!ZZZ bugfix to be done
          call wiso_flxoce(2,rbot(n),zbot(n),s16O(n),ts(n),r16O(n),ustar,re,ssq, evap_16O(n),&
               qbot(n),evap(n))
          call wiso_flxoce(3,rbot(n),zbot(n),sHDO(n),ts(n),rHDO(n),ustar,re,ssq, evap_HDO(n),&
               qbot(n),evap(n))
          call wiso_flxoce(4,rbot(n),zbot(n),s18O(n),ts(n),r18O(n),ustar,re,ssq, evap_18O(n),&
               qbot(n),evap(n))

          !------------------------------------------------------------
          ! compute diagnostics: 2m ref T & Q, 10m wind speed squared
          !------------------------------------------------------------

          if (ocn_surface_flux_scheme .eq. 0) then ! use Large algorithm

             hol = hol*ztref/zbot(n)
             xsq = max( 1.0_r8, sqrt(abs(1.0_r8-16.0_r8*hol)) )
             xqq = sqrt(xsq)
             psix2   = -5.0_r8*hol*stable + (1.0_r8-stable)*psixhu(xqq)
             fac     = (rh/shr_const_karman) * (alz + al2 - psixh + psix2 )
             tref(n) = thbot(n) - delt*fac
             tref(n) = tref(n) - 0.01_r8*ztref   ! pot temp to temp correction
             fac     = (re/shr_const_karman) * (alz + al2 - psixh + psix2 )
             qref(n) =  qbot(n) - delq*fac

             duu10n(n) = u10n*u10n ! 10m wind speed squared

          else if (ocn_surface_flux_scheme .eq. 1) then! use COARE algorithm

             tref(n) = trf
             qref(n) = qrf
             duu10n(n) = urf**2+vrf**2
             u10n = sqrt(duu10n(n))
          endif

          if (flux_diurnal) then

             !------------------------------------------------------------
             ! update new prognostic variables
             !------------------------------------------------------------

             warm  (n) = DTiter
             salt  (n) = DSiter
             speed (n) = DViter

             if (ltwopm) then
                tSkin_day(n) = tSkin(n)
                warmmax(n) = max(DTiter,0.0_r8)
             endif

             if (ltwoam) then
                tSkin_night(n) = tSkin(n)
                cSkin_night(n) = cSkin(n)
             endif

             if ((lmidnight).and.(lfullday)) then
                qSolAvg(n) = qSolInc(n)/real(nsum+1,r8)
                windAvg(n) = windInc(n)/real(nsum+1,r8)
                ! warmMax(n) = max(DTiter,warmMaxInc(n))
                windMax(n) = max(u10n,windMaxInc(n))

                nsum = 0

                qSolInc(n) = Qsol
                windInc(n) = u10n

                ! warmMaxInc(n) = 0.0_r8
                windMaxInc(n) = 0.0_r8

                !               tSkin_night(n) = tSkin(n)
                !               cSkin_night(n) = cSkin(n)

             else

                if ((lmidnight).and.(.not.(lfullday))) then

                   nsum = 0

                   qSolInc(n) = Qsol
                   windInc(n) = u10n

                   ! warmMaxInc(n) = 0.0_r8
                   windMaxInc(n) = 0.0_r8

                else

                   nsum = nsum + 1

                   ! warmMaxInc (n) = max(DTiter,warmMaxInc(n))
                   windMaxInc (n) = max(u10n, windMaxInc(n))
                   ! windMaxInc (n) = max(Qsol, windMaxInc(n))
                   qSolInc    (n) = qSolInc(n)+Qsol
                   windInc    (n) = windInc(n)+u10n

                endif
             endif

             nInc(n) = real(nsum,r8) ! set nInc to incremented or reset nsum


             if (present(ustar_sv)) ustar_sv(n) = ustar
             if (present(re_sv   )) re_sv   (n) = re
             if (present(ssq_sv  )) ssq_sv  (n) = ssq

          else              ! mask = 0

             !------------------------------------------------------------
             ! no valid data here -- out of domain
             !------------------------------------------------------------
             warm       (n) = spval ! NEW
             salt       (n) = spval ! NEW
             speed      (n) = spval ! NEW
             regime     (n) = spval ! NEW
             tBulk      (n) = spval ! NEW
             tSkin      (n) = spval ! NEW
             tSkin_night(n) = spval ! NEW
             tSkin_day  (n) = spval ! NEW
             cSkin      (n) = spval ! NEW
             cSkin_night(n) = spval ! NEW
             warmMax    (n) = spval ! NEW
             windMax    (n) = spval ! NEW
             qSolAvg    (n) = spval ! NEW
             windAvg    (n) = spval ! NEW
             warmMaxInc (n) = spval ! NEW
             windMaxInc (n) = spval ! NEW
             qSolInc    (n) = spval ! NEW
             windInc    (n) = spval ! NEW
             nInc       (n) = 0.0_r8 ! NEW

             sen   (n)    = spval  ! sensible         heat flux  (W/m^2)
             lat   (n)    = spval  ! latent           heat flux  (W/m^2)
             lwup  (n)    = spval  ! long-wave upward heat flux  (W/m^2)
             evap  (n)    = spval  ! evaporative water flux ((kg/s)/m^2)
             evap_16O (n) = spval  ! water tracer flux (kg/s)/m^2)
             evap_HDO (n) = spval  ! HDO tracer flux  (kg/s)/m^2)
             evap_18O (n) = spval  ! H218O tracer flux (kg/s)/m^2)
             taux  (n)    = spval  ! x surface stress (N)
             tauy  (n)    = spval  ! y surface stress (N)
             tref  (n)    = spval  ! 2m reference height temperature (K)
             qref  (n)    = spval  ! 2m reference height humidity (kg/kg)
             duu10n(n)    = spval  ! 10m wind speed squared (m/s)^2

             if (present(ustar_sv)) ustar_sv(n) = spval
             if (present(re_sv   )) re_sv   (n) = spval
             if (present(ssq_sv  )) ssq_sv  (n) = spval

          endif   ! mask

       endif ! flux diurnal logic

    ENDDO ! end n loop

  END subroutine flux_atmOcn_diurnal

  !===============================================================================
  ! !IROUTINE: shr_flux_MOstability -- Monin-Obukhov BL stability functions
  !
  ! !DESCRIPTION:
  !
  !    Monin-Obukhov boundary layer stability functions, two options:
  !    turbulent velocity scales or gradient and integral functions
  !    via option = shr_flux_MOwScales or shr_flux_MOfunctions
  !
  ! !REVISION HISTORY:
  !    2007-Sep-19 - B. Kauffman, Bill Large - first version
  !===============================================================================
  subroutine flux_MOstability(logunit,option,arg1,arg2,arg3,arg4,arg5)

    ! !USES:

    implicit none

    ! !INPUT/OUTPUT PARAMETERS:
    integer    ,intent(in)           :: logunit
    integer    ,intent(in)           :: option ! shr_flux_MOwScales or MOfunctions
    real(r8)   ,intent(in)           :: arg1   ! scales: uStar (in)  funct: zeta (in)
    real(r8)   ,intent(inout)        :: arg2   ! scales: zkB   (in)  funct: phim (out)
    real(r8)   ,intent(out)          :: arg3   ! scales: phim  (out) funct: phis (out)
    real(r8)   ,intent(out)          :: arg4   ! scales: phis  (out) funct: psim (out)
    real(r8)   ,intent(out),optional :: arg5   ! scales:    (unused) funct: psis (out)

    !----- local variables -----
    real(r8)           :: zeta  ! z/L
    real(r8)           :: uStar ! friction velocity
    real(r8)           :: zkB   ! (height)*(von Karman)*(surface bouyancy flux)
    real(r8)           :: phim  ! momentum    gradient function or scale
    real(r8)           :: phis  ! temperature gradient function or scale
    real(r8)           :: psim  ! momentum    integral function or scale
    real(r8)           :: psis  ! temperature integral function or scale
    real(r8)           :: temp  ! temporary-variable/partial calculation

    !----- local variables, stable case -----
    real(r8),parameter :: uStarMin = 0.001_r8 ! lower bound on uStar
    real(r8),parameter :: a = 1.000_r8  ! constant from Holtslag & de Bruin, equation 12
    real(r8),parameter :: b = 0.667_r8  ! constant from Holtslag & de Bruin, equation 12
    real(r8),parameter :: c = 5.000_r8  ! constant from Holtslag & de Bruin, equation 12
    real(r8),parameter :: d = 0.350_r8  ! constant from Holtslag & de Bruin, equation 12

    !----- local variables, unstable case -----
    real(r8),parameter :: a2 = 3.0_r8   ! constant from Wilson, equation 10

    !----- formats -----
    character(*),parameter :: subName = '(shr_flux_MOstability) '
    character(*),parameter ::   F00 = "('(shr_flux_MOstability) ',4a)"
    character(*),parameter ::   F01 = "('(shr_flux_MOstability) ',a,i5)"

    !-------------------------------------------------------------------------------
    ! Notes::
    !   o this could be two routines, but are one to help keep them aligned
    !   o the stable calculation is taken from...
    !     A.A.M. HoltSlag and H.A.R. de Bruin, 1988:
    !     "Applied Modeling of the Nighttime Surface Energy Balance over Land",
    !     Journal of Applied Meteorology, Vol. 27, No. 6, June 1988, 659-704
    !   o the unstable calculation is taken from...
    !     D. Keith Wilson, 2001: "An Alternative Function for the Wind and
    !     Temperature Gradients in Unstable Surface Layers",
    !     Boundary-Layer Meteorology, 99 (2001), 151-158
    !-------------------------------------------------------------------------------

    !----- check for consistancy between option and arguments ------------------
    if (debug > 1) then
       if (debug > 2) write(logunit,F01) "enter, option = ",option
       if ( option == shr_flux_MOwScales .and. present(arg5) ) then
          write(logunit,F01) "ERROR: option1 must have four arguments"
          call shr_sys_abort(subName//"option inconsistant with arguments")
       else if ( option == shr_flux_MOfunctions .and. .not. present(arg5) ) then
          write(logunit,F01) "ERROR: option2 must have five arguments"
          call shr_sys_abort(subName//"option inconsistant with arguments")
       else
          write(logunit,F01) "invalid option = ",option
          call shr_sys_abort(subName//"invalid option")
       end if
    end if

    !------ velocity scales option ----------------------------------------------
    if (option == shr_flux_MOwScales) then

       !--- input ---
       uStar = arg1
       zkB   = arg2

       if (zkB >= 0.0_r8) then ! ----- stable -----
          zeta = zkB/(max(uStar,uStarMin)**3)
          temp = exp(-d*zeta)
          phim = uStar/(1.0_r8 + zeta*(a + b*(1.0_r8 + c - d*zeta)*temp))
          phis = phim
       else                    ! ----- unstable -----
          temp = (zkB*zkB)**(1.0_r8/a2)   ! note: zkB < 0, zkB*zkB > 0
          phim = sqrt(uStar**2 + shr_flux_MOgammaM*temp)
          phis = sqrt(uStar**2 + shr_flux_MOgammaS*temp)
       end if

       !--- output ---
       arg3 = phim
       arg4 = phis
       !  arg5 = <unused>

       !------ stability function option -------------------------------------------
    else if (option == shr_flux_MOfunctions) then

       !--- input ---
       zeta  = arg1

       if (zeta >= 0.0_r8) then ! ----- stable -----
          temp = exp(-d*zeta)
          phim =        1.0_r8 + zeta*(a + b*(1.0_r8 + c - d*zeta)*temp)
          phis = phim
          psim = -a*zeta - b*(zeta - c/d)*temp - b*c/d
          psis = psim
       else                    ! ----- unstable ----
          temp = (zeta*zeta)**(1.0_r8/a2)   ! note: zeta < 0, zeta*zeta > 0
          phim = 1.0_r8/sqrt(1.0_r8 + shr_flux_MOgammaM*temp)
          phis = 1.0_r8/sqrt(1.0_r8 + shr_flux_MOgammaS*temp)
          psim = a2*log(0.5_r8 + 0.5_r8/phim)
          psis = a2*log(0.5_r8 + 0.5_r8/phis)
       end if

       !--- output ---
       arg2 = phim
       arg3 = phis
       arg4 = psim
       arg5 = psis
       !----------------------------------------------------------------------------
    else
       write(logunit,F01) "invalid option = ",option
       call shr_sys_abort(subName//"invalid option")
    endif

  end subroutine flux_MOstability

  !===============================================================================
  ! !DESCRIPTION:
  !
  !   COARE v3.0 parametrisation
  !
  ! !REVISION HISTORY:
  !   2013-Nov-22: Thomas Toniazzo's adaptation of Chris Fairall's code,
  !    downloaded from
  !    ftp://ftp1.esrl.noaa.gov/users/cfairall/wcrp_wgsf/computer_programs/cor3_0/
  !     * no wave, standard coare 2.6 charnock
  !     * skin parametrisation also off (would require radiative fluxes and
  !      rainrate in input)
  !     * added diagnostics, comments and references
  !===============================================================================
  subroutine cor30a(ubt,vbt,tbt,qbt,rbt        &    ! in atm params
       & ,uss,vss,tss,qss            &    ! in surf params
       & ,zbl,zbu,zbt,zrfu,zrfq,zrft &    ! in heights
       & ,tau,hsb,hlb                &    ! out: fluxes
       & ,zo,zot,zoq,L,usr,tsr,qsr   &    ! out: ss scales
       & ,Cd,Ch,Ce                   &    ! out: exch. coeffs
       & ,trf,qrf,urf,vrf)                ! out: reference-height params

    ! !USES:

    IMPLICIT NONE

    ! !INPUT/OUTPUT PARAMETERS:

    real(r8),intent(in) :: ubt,vbt,tbt,qbt,rbt,uss,vss,tss,qss
    real(r8),intent(in) :: zbl,zbu,zbt,zrfu,zrfq,zrft
    real(r8),intent(out):: tau,hsb,hlb,zo,zot,zoq,L,usr,tsr,qsr,Cd,Ch,Ce &
         & ,trf,qrf,urf,vrf

    real(r8) ua,va,ta,q,rb,us,vs,ts,qs,zi,zu,zt,zq,zru,zrq,zrt ! internal vars

    real(r8):: cpa,rgas,grav,pi,von,beta ! phys. params
    real(r8):: le,rhoa,cpv               ! derived phys. params
    real(r8):: t,visa,du,dq,dt           ! params of problem

    real(r8):: u10,zo10,zot10,cd10,ch10,ct10,ct,cc,ribu,zetu,l10,charn ! init vars
    real(r8):: zet,rr,bf,ug,ut     ! loop iter vars
    real(r8):: cdn_10,chn_10,cen_10  ! aux. output vars

    integer(in) i, nits ! iter loop counters

    integer(in) jcool                  ! aux. cool-skin vars
    real(r8) dter, wetc, dqer

    ua=ubt  !wind components (m/s) at height zu (m)
    va=vbt
    ta=tbt  !bulk air temperature (K), height zt
    Q =qbt  !bulk air spec hum (kg/kg), height zq
    rb=rbt  ! air density
    us=uss  !surface current components (m/s)
    vs=vss
    ts=tss  !bulk water temperature (K) if jcool=1, interface water T if jcool=0
    qs=qss  !bulk water spec hum (kg/kg) if jcool=1 etc
    zi=zbl  !PBL depth (m)
    zu=zbu  !wind speed measurement height (m)
    zt=zbt  !air T measurement height (m)
    zq=zbt  !air q measurement height (m)
    zru=zrfu ! reference height for st.diagn.U
    zrq=zrfq ! reference height for st.diagn.T,q
    zrt=zrft ! reference height for st.diagn.T,q

    !**** constants
    Beta= 1.2_r8
    von = 0.4_r8
    pi  = 3.141593_r8
    grav= SHR_CONST_G
    Rgas= SHR_CONST_RGAS
    cpa = SHR_CONST_CPDAIR

    !*** physical parameters
    Le  = SHR_CONST_LATVAP -.00237e6_r8*(ts-273.16_r8)
    !   cpv = shr_const_cpdair*(1.0_r8 + shr_const_cpvir*Qs) ! form in NCAR code
    cpv = cpa*(1.0_r8+0.84_r8*Q)
    !   rhoa= P/(Rgas*ta*(1+0.61*Q)) ! if input were pressure
    rhoa= rb

    ! parametrisation for air kinematic viscosity (Andreas 1989,p.31)
    t   = ta-273.16_r8
    visa= 1.326e-5_r8*(1.0_r8+6.542e-3_r8*t+8.301e-6_r8*t*t-4.84e-9_r8*t*t*t)

    du  = sqrt((ua-us)**2+(va-vs)**2)
    dt  = ts-ta -.0098_r8*zt
    dq  = Qs-Q

    !*** don't use cool-skin params for now, but assign values to Ter and Qer
    jcool=0_IN
    dter=0.3_r8
    wetc=0.622_r8*Le*Qs/(Rgas*ts**2)
    dqer=wetc*dter

    !***************** Begin bulk-model calculations ***************

    !*************** first guess
    ug=0.5_r8

    ut   = sqrt(du*du+ug*ug)
    u10  = ut*log(10.0_r8/1.0e-4_r8)/log(zu/1.0e-4_r8)
    usr  = .035_r8*u10
    zo10 = 0.011_r8*usr*usr/grav+0.11_r8*visa/usr
    Cd10 = (von/log(10.0_r8/zo10))**2
    Ch10 = 0.00115_r8
    Ct10 = Ch10/sqrt(Cd10)
    zot10= 10.0_r8/exp(von/Ct10)
    Cd   =(von/log(zu/zo10))**2
    Ct   = von/log(zt/zot10)
    CC   = von*Ct/Cd

    ! Bulk Richardson number
    Ribu=-grav*zu/ta*((dt-dter*jcool)+.61_r8*ta*dq)/ut**2
    ! initial guess for stability parameter...
    if (Ribu .LT. 0.0_r8) then
       ! pbl-height dependent
       zetu=CC*Ribu/( 1.0_r8 - (.004_r8*Beta**3*zi/zu) * Ribu )
    else
       zetu=CC*Ribu*(1.0_r8 + 27.0_r8/9.0_r8*Ribu/CC)
    endif
    ! ...and MO length
    L10=zu/zetu

    if (zetu .GT. 50.0_r8) then
       nits=1_IN
    else
       nits=3_IN
    endif

    usr =  ut*von/(log(zu/zo10)-psiuo(zu/L10))
    tsr = (dt-dter*jcool)*von/(log(zt/zot10)-psit_30(zt/L10))
    qsr = (dq-dqer*jcool)*von/(log(zq/zot10)-psit_30(zq/L10))

    ! parametrisation for Charney parameter (section 3c of Fairall et al. 2003)
    charn=0.011_r8
    if (ut .GT. 10.0_r8) then
       charn=0.011_r8+(ut-10.0_r8)/(18.0_r8-10.0_r8)*(0.018_r8-0.011_r8)
    endif
    if (ut .GT. 18.0_r8) then
       charn=0.018_r8
    endif

    !***************  iteration loop ************
    do i=1, nits

       ! stability parameter
       zet=-von*grav*zu/ta*(tsr*(1.0_r8+0.61_r8*Q)+.61_r8*ta*qsr)/(usr*usr)/(1.0_r8+0.61_r8*Q)

       ! momentum roughness length...
       zo = charn*usr*usr/grav+0.11_r8*visa/usr
       ! ...& MO length
       L  = zu/zet

       ! tracer roughness length
       rr = zo*usr/visa
       zoq= min(1.15e-4_r8,5.5e-5_r8/rr**.6_r8)
       zot= zoq ! N.B. same for vapour and heat

       ! new surface-layer scales
       usr =  ut            *von/(log(zu/zo )-psiuo(zu/L))
       tsr = (dt-dter*jcool)*von/(log(zt/zot)-psit_30(zt/L))
       qsr = (dq-dqer*jcool)*von/(log(zq/zoq)-psit_30(zq/L))

       ! gustiness parametrisation
       Bf=-grav/ta*usr*(tsr+.61_r8*ta*qsr)
       if (Bf .GT. 0.0_r8) then
          ug=Beta*(Bf*zi)**.333_r8
       else
          ug=.2_r8
       endif
       ut=sqrt(du*du+ug*ug)

    enddo
    !***************     end loop    ************

    !******** fluxes @ measurement heights zu,zt,zq ********
    tau= rhoa*usr*usr*du/ut                !stress magnitude
    hsb=-rhoa*cpa*usr*tsr                  !heat downwards
    hlb=-rhoa*Le*usr*qsr                   !wv downwards

    !****** transfer coeffs relative to ut @meas. hts ******
    Cd= tau/rhoa/ut/max(.1_r8,du)
    if (tsr.ne.0._r8) then
       Ch= usr/ut*tsr/(dt-dter*jcool)
    else
       Ch= usr/ut* von/(log(zt/zot)-psit_30(zt/L))
    endif
    if (qsr.ne.0.0_r8) then
       Ce= usr/ut*qsr/(dq-dqer*jcool)
    else
       Ce= usr/ut* von/(log(zq/zoq)-psit_30(zq/L))
    endif

    !**********  10-m neutral coeff relative to ut *********
    Cdn_10=von*von/log(10.0_r8/zo)/log(10.0_r8/zo)
    Chn_10=von*von/log(10.0_r8/zo)/log(10.0_r8/zot)
    Cen_10=von*von/log(10.0_r8/zo)/log(10.0_r8/zoq)

    !**********  reference-height values for u,q,T *********
    urf=us+(ua-us)*(log(zru/zo)-psiuo(zru/L))/(log(zu/zo)-psiuo(zu/L))
    vrf=vs+(va-vs)*(log(zru/zo)-psiuo(zru/L))/(log(zu/zo)-psiuo(zu/L))
    qrf=qs-dq*(log(zrq/zoq)-psit_30(zrq/L))/(log(zq/zoq)-psit_30(zq/L))
    trf=ts-dt*(log(zrt/zot)-psit_30(zrt/L))/(log(zt/zot)-psit_30(zt/L))
    trf=trf+.0098_r8*zrt

  end subroutine cor30a

  !===============================================================================
  ! !IROUTINE: PSIUo
  !
  ! !DESCRIPTION:
  !
  !   momentum stability functions adopted in COARE v3.0 parametrisation.
  !   Chris Fairall's code (see cor30a)
  !
  ! !REVISION HISTORY:
  !   22/11/2013: Thomas Toniazzo: comments added
  !===============================================================================

  real (r8) function psiuo(zet)
    ! !INPUT/OUTPUT PARAMETERS:
    real(r8),intent(in)  :: zet
    real(r8) ::c,x,psik,psic,f
    !-----------------------------------------------------------------
    ! N.B.: z0/L always neglected compared to z/L and to 1
    !-----------------------------------------------------------------
    if(zet>0.0_r8)then
       ! Beljaars & Holtslag (1991)
       c=min(50._r8,.35_r8*zet)
       psiuo=-((1.0_r8+1.0_r8*zet)**1.0_r8+.667_r8*(zet-14.28_r8)/exp(c)+8.525_r8)
    else
       ! Dyer & Hicks (1974) for weak instability
       x=(1.0_r8-15.0_r8*zet)**.25_r8                   ! 15 instead of 16
       psik=2.0_r8*log((1.0_r8+x)/2.0_r8)+log((1.0_r8+x*x)/2.0_r8)-2.0_r8*atan(x)+2.0_r8*atan(1.0_r8)
       ! Fairall et al. (1996) for strong instability (Eq.(13))
       x=(1.0_r8-10.15_r8*zet)**.3333_r8
       psic= 1.5_r8*log((1.0_r8+x+x*x)/3.0_r8)-sqrt(3.0_r8)*atan((1.0_r8+2.0_r8*x)/sqrt(3.0_r8)) &
            & +4.0_r8*atan(1.0_r8)/sqrt(3.0_r8)
       f=zet*zet/(1.0_r8+zet*zet)
       psiuo=(1.0_r8-f)*psik+f*psic
    endif
  END FUNCTION psiuo

  !===============================================================================
  ! !IROUTINE: PSIT_30
  !
  ! !DESCRIPTION:
  !
  !   momentum stability functions adopted in COARE v3.0 parametrisation.
  !   Chris Fairall's code (see cor30a)
  !
  ! !REVISION HISTORY:
  !   22/11/2013: Thomas Toniazzo: comments added
  !===============================================================================
  real (r8) function psit_30(zet)
    ! !INPUT/OUTPUT PARAMETERS:
    real(r8),intent(in)  :: zet
    ! !EOP
    real(r8) ::c,x,psik,psic,f
    !-----------------------------------------------------------------
    ! N.B.: z0/L always neglected compared to z/L and to 1
    !-----------------------------------------------------------------
    if(zet>0.0_r8)then
       ! Beljaars & Holtslag (1991)
       c=min(50._r8,.35_r8*zet)
       psit_30=-((1.0_r8+2.0_r8/3.0_r8*zet)**1.5_r8+.667_r8*(zet-14.28_r8)/exp(c)+8.525_r8)
    else
       ! Dyer & Hicks (1974) for weak instability
       x=(1.0_r8-15.0_r8*zet)**.5_r8                    ! 15 instead of 16
       psik=2.0_r8*log((1.0_r8+x)/2.0_r8)
       ! Fairall et al. (1996) for strong instability
       x=(1.0_r8-(34.15_r8*zet))**.3333_r8
       psic= 1.5_r8*log((1.0_r8+x+x*x)/3.0_r8)-sqrt(3.0_r8)*atan((1.0_r8+2.0_r8*x)/sqrt(3.0_r8)) &
            & +4.0_r8*atan(1.0_r8)/sqrt(3.0_r8)
       f=zet*zet/(1.0_r8+zet*zet)
       psit_30=(1.0_r8-f)*psik+f*psic
    endif
  end FUNCTION psit_30

end module shr_flux_mod
